5. Численное дифференцирование Задача приближенного вычисления производной может возникнуть в тех случаях, когда неизвестно аналитическое выражение для исследуемой функции. Функция может быть задана таблично, или известен только график функции, полученный, например, в результате показаний датчиков параметров технологического процесса. Иногда, при решении некоторых задач на компьютере, из-за громоздкости выкладок может оказаться более удобным вычисление производных численным методом, чем аналитическим. При этом, разумеется, необходимо обосновать применяемый численный метод, т.е. убедиться в том, что погрешность численного метода находится в приемлемых границах. Одним из эффективных методов решения дифференциальных уравнений является разностный метод, когда вместо искомой функции рассматривается таблица её значений в определенных точках, при этом производные приближенно заменяются разностными формулами

5.1.Графическое дифференцирование Пусть известен график функции y = f(x) на отрезке [a, b]. Можно построить график производной функции, вспомнив ее геометрический смысл. Воспользуемся тем фактом, что производная функции в точке x равна тангенсу угла наклона к оси абсцисс касательной к её графику в этой точке абсцисс. Если x = x 0, то найдем y 0 = f(x 0 ) с помощью графика и затем проведем касательную AB к графику функции в точке (x 0, y 0 ) (рис. 5.1). Проведем прямую, параллельную касательной AB, через точку (–1, 0) и найдем точку y 1 её пересечения с осью ординат. Тогда значение y 1 равно тангенсу угла наклона касательной к оси абсцисс, т.е. производной функции f(x) в точке x 0 : и точка M 0 (x 0, y 1 ) принадлежит графику производной. Рис. 5.1 x0x0 –10 y0y0 y x y1y1 M0M0 A B α

Чтобы построить график производной необходимо разбить отрезок [a, b] на несколько частей точками x i, затем для каждой точки графически построить значение производной и соединить полученные точки плавной кривой с помощью лекал. На рис. 5.2 показано построение пяти точек M 1, M 2, …, M 5 и графика производной. Алгоритм построения графика производной: 1. Строим касательную к графику функции y = f(x) в точке (x 1, f(x 1 )); из точки (–1, 0) параллельно касательной в точке (x 1, f(x 1 )) проведем прямую до пересечения с осью ординат; эта точка пересечения дает значение производной. Строим точку. 2. Аналогично построим остальные точки M 2, M 3, M 4 и M Соединяем точки M 1, M 2, M 3, M 4, M 5 плавной кривой. Полученная кривая является графиком производной.

Рис. 5.2 Точность графического способа определения производной невысока. Мы приводим описание этого способа только в учебных целях. Замечание. Если в алгоритме построения графика производной вместо точки (–1, 0) взять точку (–l, 0), где l > 0, то график будет построен в другом масштабе по оси ординат. x1x1 –10 y x M4M4 x2x2 M1M1 x4x4 M2M2 x3x3 M3M3 x5x5 M5M5

5.2. Разностные формулы Разностные формулы для приближенного вычисления производной подсказаны самим определением производной. Пусть значения функции в точках x i обозначены через y i : y i = f(x i ), x i = a + ih, i = 0, 1, …, n; h = (b – a)/n. Мы рассматриваем случай равномерного распределения точек на отрезке [a, b]. Для приближенного вычисления производных в точках x i можно использовать следующие разностные формулы, или разностные производные. (5.1) (5.2) (5.3) Разностные формулы для обыкновенных производных

Так как предел отношения (5.1) при h 0 равен правой производной в точке x i, то это отношение иногда называют правой разностной производной в точке x i. По аналогичной причине отношение (5.2) называют левой разностной производной в точке x i. Отношение (5.3) называют центральной разностной производной в точке x i. Оценим погрешность разностных формул (5.1) (5.3), предполагая, что функция f(x) разлагается в ряд Тейлора в окрестности точки x i : (5.4) Полагая в (5.4) x = x i + h или x = x i – h, получим (5.5) (5.6) Учитывая (5.5) и (5.6) имеем (5.7) (5.8)

(5.9) Из последних соотношений следует, что разностная формула (5.3) имеет погрешность на порядок меньшую, чем разностные формулы (5.1) и (5.2). Пример 5.1. Вычислить приближенно производные с помощью разностных формул и сравнить с точными значениями производной функции y = sin x в точках отрезка [0, 1]. Решение в Excel. Составим таблицу значений функции y = sin πx на отрезке [0, 1] с шагом 0,1 (табл. 5.1). В ячейку B2 введем формулу =SIN(A2*3, ) и протянем ячейку B2 маркером заполнения до ячейки B12. В ячейки C3, D3, E3 и F3 введем соответственно формулы =(B4-B3)/(A4-A3), =(B3-B2)/(A3-A2), =(B4-B2)/(A4-A2) и =3, *COS(A3*3, ), выделим диапазон C3:F3 и протянем маркером заполнения до строки 11. Результат представлен также графиком на рис Как видно на рисунке, разностная формула (5.3) дает практически те же значения, что и формула точной производной.

Табл. 5.1 Производные высших порядков можно приближенно вычислять по формулам, полученным с помощью последовательного применения разностных соотношений (5.1) (5.3). Разностная формула для второй производной (разностная производная второго порядка) имеет вид: (5.10) Непосредственной подстановкой разложений (5.5) и (5.6) в формулу (5.10) можно получить зависимость между второй производной функции и разностной формулой для производной второго порядка: (5.11)

Рис. 5.3.

Разностные формулы для частных производных Разностные формулы для частных производных аналогичны формулам (5.1) (5.3), (5.10). Пусть функция двух переменных определена в прямоугольной области a 1 x b 1, a 2 y b 2. Определение 5.4. Назовем сеточной областью множество точек (x i, y j ), где x i = a 1 + ih 1, i = 0, 1, …, n 1 ; h 1 = (b 1 – a 1 )/n 1, y j = a 2 + jh 2, j = 0, 1, …, n 2 ; h 2 = (b 2 – a 2 )/n 2. На рис. 5.3 изображена сеточная область для n 1 = 5, n 2 = 4. Эта сеточная область состоит из 30 точек. Рис. 5.3 a1a1 y x a2a2 b2b2 b1b1

Введем обозначение u i, j = f(x i, y j ). Тогда для частных производных первого и второго порядка по переменной x можно записать разностные формулы (5.1) (5.3) и (5.10): (5.12) (5.13) (5.14) (5.15) Аналогичные разностные формулы можно записать и для частных производных первого и второго порядка по переменной y: (5.16) (5.17)

(5.18) (5.19) Запишем разностные формулы для смешанных производных: (5.20) (5.21) (5.22)

5.3.Вычисление производных с помощью интерполяционных формул Разностные формулы для производных не годятся в тех случаях, когда необходимо вычислить производную в произвольной точке, не совпадающей с узлами x i. В этих случаях естественно воспользоваться интерполяционными формулами Равномерное распределение узлов Пусть заданы значения функции y i = f(x i ) в узлах x i, равномерно распределенных на отрезке [a, b]: x i = a + ih, i = 0, 1, …, n; h = (b – a)/n. Построим формулы приближенного дифференцирования с помощью первой интерполяционной формулы Ньютона (4.13), которую после некоторых упрощений можно записать в виде (5.23) где q = (x – x 0 )/h. Мы здесь предполагаем, что узлы интерполяции распределены на отрезке [a, b] равномерно с шагом h. Выберем в качестве x 0 узел, ближайший к точке x. Дифференцируя (5.23) по переменной x и учитывая, что

получим приближенные формулы для вычисления производных: (5.24) (5.25) Чем больше слагаемых в формулах (5.24) и (5.25) используется, тем выше точность вычисления производных. Эти формулы дают хорошие приближения для точек, близких к значению x 0 (левому концу отрезка [a, b]). Если x = x 0, то q = 0 и для вычисления производных в узлах x i получим формулы (5.26) (5.27) Очевидно, что вторая интерполяционная формула Ньютона позволяет вывести формулы для вычисления производных в точках, близких к правому краю x n отрезка [a, b]:

(5.28) (5.29) (5.30) При x = x n получим (5.31) (5.32)

Неравномерное распределение узлов Пусть известны значения y i функции в узлах x i : x 0 < x 1 < x 2 < … < x n – 1 < x n. В этом случае для вычисления производных используют интерполяционный многочлен. Предположим, что точка x расположена ближе к начальному узлу x 0, тогда мы можем применить первый интерполяционный многочлен Ньютона Здесь использованы обозначения для разделенных разностей:

Обозначим через ξ i разность (x – x i ) и запишем многочлен Ньютона в виде Теперь можем вывести формулы для производных: (5.38) Оставляя в (5.38) несколько слагаемых, получим формулы для приближенного вычисления производных. При этом порядок погрешности по отношению к шагу разбиения h равен числу оставленных членов, или разности между числом узлов интерполяции и порядком производной. Приведем некоторые простые формулы (h = max h i ): 1. Первая производная по двум точкам: (5.39)

Первая производная по трем точкам: (5.40) 2. Вторая производная по трем точкам: (5.41) Вторая производная по четырем точкам: (5.42)

3. Третья производная по четырем точкам: (5.43) Пример 5.5. Вычислить производную второго порядка функции y = sin πx в точке x = 0,4, пользуясь таблицей значений (см. табл.5.8). Решение в Excel. Воспользуемся формулами (5.41) и (5.42). Для формулы (5.41) необходимы только три точки. Создадим в программе Excel макрос функцию для вычисления производной по формуле (5.41). Для этого выберем меню «Сервис Макрос Редактор Visual Basic» и в открывшемся окне выполним команду «Insert Module», затем введем: Function pr_541(x0, x1, x2, y0, y1, y2) pr_541 = 2 * ((y0 - y1) / (x0 - x1) - (y1 - y2) / (x1 - x2)) / (x0 - x2) End Function Перейдем в Excel и в ячейке D4 запишем формулу =pr_541(B2;B3;B4;C2;C3;C4), и скопируем D4 в D5, D6. Таблица 5.8 Наиболее точно вторую производную приближает значение в D6. Это объясняется тем, что в формуле, записанной в D6 «=pr_541(B4;B5;B6;C4;C5;C6)» точка x = 0,4 близка к начальному узлу x0 =0,3.

5.4.Практическая оценка погрешности. Метод Рунге-Ромберга Рассмотрим метод Рунге для повышения порядка точности формул для вычисления производных на примере формулы (5.11): (5.44) Предполагая, что функция четырежды дифференцируема, погрешность разностной формулы для второй производной можно представить в виде (5.45) Вычислим по формуле (5.44) вторую производную в той же точке x i, но используя сетку с шагом 2h: (5.46) Погрешность формулы (5.46) имеет вид: (5.47) Вычитая (5.45) из (5.47), получим (5.48)

Отсюда для погрешности получим первую формулу Рунге: (5.49) Теперь вторую производную в точке x i можно вычислить по уточненной формуле (второй формуле Рунге): (5.50)