Простейшие формулы численного дифференцирования

Лекция №3

Численное дифференцирование

Интерполяционный полином Ньютона

Методы численного дифференцирования широко используются на практике, когда численно необходимо решить одно или несколько дифференциальных уравнений. Методы численного дифференцирования привлекаются в тех случаях, когда необходимо найти ту или иную производные, поиск которой аналитически затруднен или просто невозможно, если функция задана таблично. В этих случаях исходную функцию y = f(x) заменяют, аппроксимируют более простой функцией f(x) и считают, что приближенно верно равенство: y¢(x) » f¢(x). При построении аппроксимирующей функции могут быть использованы разнообразные методы, в том числе и те, которые представлены в предыдущей лекции №2.

Рассмотрим еще один способ аппроксимации, который не рассмотрен в предыдущее лекции, — аппроксимация интерполяционным многочленом Ньютона.

Пусть функция y = f(x) определена таблично, т.е. задан набор пар (x₁,y₁), …, (x_n,y_n) такой, что y_i  = f(x_i), i = 1,…,n. Определим так называемые разделенные разности функции y = f(x):

f(x_i,x_j) = [f(x_i) - f(x_j)]/(x_i - x_j),

f(x_i,x_j,x_k) = [f(x_i,x_j) - f(x_j,x_k)]/(x_i - x_k),                                             (1)

f(x_i,x_j,x_k,x_l) = [f(x_i,x_j,x_k) - f(x_j,x_k,x_l)]/(x_i - x_l)

и т.д. В (1) определены разделенные разности первого, второго и третьего порядков.

Пусть P(x) — многочлен степени n -1 изучим его разделенные разности. Рассмотрим разность P(x) - P(x₁), она обращается в нуль при x = x₁ и поэтому делится на x - x₁. Таким образом, первая разделенная разность P(x,x₁) = [P(x) - P(x₁)]/(x - x₁) многочлена степени n -1 является полиномом степени n -2. Аналогично вторая разделенная разность P(x,x₁,x₂) является полиномом степени n -3 и т.д. Этот процесс рассуждений можно продолжить вплоть до разделенной разности P(x,x₁,…,x_{n -1}), которая является полиномом нулевой степени, т.е. константой. Более высокие разделенные разности равны тождественно нулю.

Перепишем (1) применительно к полиному P(x) выбрав в качестве первого аргумента x, тогда

P(x) = P(x₁) + (x - x₁)P(x,x₁),

P(x,x₁) = P(x₁,x₂) + (x - x₂)P(x,x₁,x₂),                                           (2)

P(x,x₁,x₂) = P(x₁,x₂,x₃) + (x - x₃)P(x,x₁,x₂,x₃)

и т.д. Подставляя выражения (2) друг в друга, получим формулу

(3)

по которой многочлен степени n -1 выражается через свои разделенные разности значений в узлах x₁,…,x_n.

Подставляя теперь в (3) разделенные разности аппроксимируемой функции, получим интерполяционную формулу Ньютона

            (4)

Для оценки погрешности R(x) = f(x) - f(x) интерполяционной формулы Ньютона (4) можно воспользоваться тем же методом, который был использован в предыдущей лекции при оценке погрешности интерполяции полиномом Лагранжа. Повторяя дословно те же рассуждения, получаем

                                                                       (5)

где w(x) = (x - x₁)…(x - x_n),  и супремум берется по интервалу между минимальным и максимальным значениями из набора x,x₁,…,x_n.

Вводя обозначения, x _i = x - x_i, i = 1,…n продифференцируем многочлен (4) несколько раз, тогда

Общая формула имеет вид:

                     (6)

Обрывая ряд, на каком либо члене получим приближенную формулу для вычисления той или иной производной. Наиболее простые и наиболее употребительные выражения получаются, когда ограничиваются первым членом, т.е.

                      (7)

Методом индукции можно также доказать следующую формулу для оценки разделенной разности n-го порядка:

                                           (8)

Для оценки погрешности формулы дифференцирования (6) будем считать шаг сетки h ( ) достаточно малым. Пусть мы используем узлы x₁,…,x_n, в этом случае погрешность  близка к первому отброшенному члену ряда, он содержит разделенную разность f(x₁,…,x_n,x_{n +1}), которая, согласно (7), (8), примерно равна f⁽ⁿ⁾(x)/n!. Итак, имеем

.                                              (9)

Сумма в (9) имеет  слагаемых так, что верна следующая цепочка оценок

                                (10)

где использована формула Стирлинга .

Согласно формуле (10), порядок точности формулы (6) по отношению к шагу сетки равен числу узлов интерполяции минус порядок производной. Минимальное число узлов интерполяции для аппроксимации k-й производной с первым порядком точности находится из уравнения: n - k = 1, т.е. равно k + 1 или на единицу больше порядка производной.

Простейшие формулы численного дифференцирования

Часто используются равномерные сетки, на которых вид формул (6), (7) заметно упрощается и точность аппроксимации производных может повыситься.

Выберем два узла x₁ < x₂ и найдем первую производную в узле x₂, т.е.

                                                        (11)

Выберем три узла x₁ < x₂ < x₃ и найдем первую и вторую производные согласно (7) в средней точке, т.е.

                                                    (12)

                                          (13)

где считается, что .

Аналогично (11) — (13) можно вывести формулы для аппроксимации производных с более высоким порядком точности.

Так, согласно (6), для первой производной в середине интервала по четырем узлам сетки имеем

              (14)

а для второй производной в центре по пяти узлам

            (15)

Обычно для априорной оценки точности формул, подобных (11) — (15), используют разложение в ряд Тейлора. Пусть функция y = f(x) имеет непрерывную четвертую производную, тогда верно следующее разложение

                         (16)

где , . Подставляя разложение (16), например, в (13), находим

            (17)

Метод разложения (16), (17) наиболее часто используется для априорной оценки порядка аппроксимации разностных схем.

 

Рис.1. Абсолютные ошибки формул численного
дифференцирования (18) — (21)

 

Проведем теперь вычислительный эксперимент по сравнению различных схем численного дифференцирования (11) — (15).

В начале сравним разные методы численной оценки первой производной функции на примере известной функции y = sin(x) на отрезке [0,p]. Перепишем (11), (12) и (14) в более удобном для нас виде с учетом того, что задана равномерная сетка 0 = x₁ < …< x_n = p., где h = x_{i +1} - x_i — шаг сетки.

Пусть вычисляется первая производная с помощью “правой” конечной разности:

;                                             (18)

с помощью “левой” конечной разности:

;                                                 (19)

с помощью центральной разности второго порядка точности:

;                                      (20)

и, наконец, с помощью формулы (14) повышенной точности — четвертого порядка:

.                     (21)

Вычисления всех четырех производных производится в программе, код которой представлен на листинге_№1. Там же подсчитываются соответствующие абсолютные погрешности и выводятся на едином графике, представленном на рис.1.

 

Листинг_№1

%Программа сравнения простейших формул

%численного дифференцирования с точным

%значением производной

%определим шаг сетки

h=0.2;

%определим сетку

x=0:h:pi;

n=length(x);

%анализируемая функция y=sin(x)

%точное значение производной dy=cos(x)

dy=cos(x);

%определим производную по формуле

%правой конечной разности и соответствующую

%абсолютную погрешность

for i=1:(n-1)

dy1(i)=(sin(x(i+1))-sin(x(i)))/h;

er1(i)=abs(dy(i)-dy1(i));

end

%определим производную по формуле

%левой конечной разности и соответствующую

%абсолютную погрешность

for i=2:n

dy2(i)=(sin(x(i))-sin(x(i-1)))/h;

er2(i)=abs(dy(i)-dy2(i));

end

%определим производную по формуле

%центральной разности и соответствующую

%абсолютную погрешность

for i=2:(n-1)

dy3(i)=(sin(x(i+1))-sin(x(i-1)))/(2*h);

er3(i)=abs(dy(i)-dy3(i));

end

%определим производную по формуле (21)

%четвертого порядка точности и соответствующую

%абсолютную погрешность, которая вычисляется

%в точке x(i)-0.5*h

for i=3:(n-1)

dy4(i)=(-sin(x(i+1))+27*sin(x(i))-...

     27*sin(x(i-1))+sin(x(i-2)))/(24*h);

er4(i)=abs(cos(x(i)-0.5*h)-dy4(i));

end

%рисуем все три абсолютные ошибки на одном графике

plot(x([1:(n-1)]),er1([1:(n-1)]),'-o',...

        x([2:n]),er2([2:n]),'-p',...

x([2:(n-1)]),er3([2:(n-1)]),'-h',...

x([3:(n-1)]),er4([3:(n-1)]),'-*');

 

Анализ рис.1 показывает значительное отличие ошибок аппроксимации первых двух схем первого порядка точности от двух других — второго и четвертого порядка точности соответственно.

Аналогично изучим путем сравнения две формулы для численной оценки второй производной функции на равномерной сетке. Перепишем формулы (13), (15) в тех обозначениях, которые уже были использованы для численной аппроксимации первой производной, т.е.

,                                                    (22)

.                 (23)

Согласно (13), (15), формула (22) аппроксимирует вторую производную со вторым порядком точности, а формула (23) — с четвертым. На листинге_№2 приведен код соответствующей программы, а на рис.2 итог работы программы. На рис.2 отчетливо видно сколь большое различие в точности аппроксимации второй производной между двумя схемами (22) и (23) имеет место.

 

Листинг_№2

%Программа двух формул для численной

%оценки второй производной

%определим шаг сетки

h=0.2;

%определим сетку

x=0:h:pi;

n=length(x);

%анализируемая функция y=sin(x)

%точное значение второй производной d2y=-sin(x)

d2y=-sin(x);

%определим вторую производную по формуле (22)

for i=2:(n-1)

d2y1(i)=(sin(x(i+1))-2*sin(x(i))+sin(x(i-1)))/h^2;

er1(i)=abs(d2y(i)-d2y1(i));

end

%определим вторую производную по формуле (23)

for i=3:(n-2)

d2y2(i)=(-sin(x(i+2))+16*sin(x(i+1))-30*sin(x(i))+...

               16*sin(x(i-1))-sin(x(i-2)))/(12*h^2);

er2(i)=abs(d2y(i)-d2y2(i));

end

%рисуем абсолютные ошибки на одном графике

plot(x([2:(n-1)]),er1([2:(n-1)]),'-o',...

x([3:(n-2)]),er2([3:(n-2)]),'-p');

 

Рис.2. Абсолютные ошибки формул численного
дифференцирования (22), (23)

 

Согласно (7) — (10), с точки зрения зависимости от шага равномерной сетки h вычислительную формулу для k-й производной можно представить в виде

                                             (24)

В (24) при вычислении z имеют место неустранимые ошибки d z, например, ошибки экспериментальных данных либо ошибки округления. Если ошибки эксперимента до известной степени допускают коррекцию с помощью мер сглаживания (например, методом наименьших квадратов), то ошибки округления являются действительно неустранимыми, т.к. разрядность чисел на компьютере ограничена. В итоге ошибку аппроксимации можно записать в виде

                                                                           (25)

из которой следует, что шаг сетки не может быть выбран как угодно малым, т.е. можно говорить об оптимальном значении шага сетки h_*. Из формулы (25) легко найти оптимальное значение шага, приравнивая нулю первую производную по h от выражения (25), т.е.

.                                                                      (26)

Дальнейшее уточнение оптимального значения шага сетки (26) затруднительно, поэтому ниже такой шаг будет обнаружен путем вычислительного эксперимента.

В конечном счете, наличие оптимального шага связано с некорректностью операции дифференцирования. Пусть функция f возмущена величиной , которая в норме ||…||_C может быть сделана как угодно малой при m ® ¥. После применения операции дифференцирования к d f, получим , и в норме ||…||_C возмущение d f¢ может быть как угодно большим при m ® ¥.

На листинге_№3 приведен код программы численной оценки второй производной функции y = sin(x) в точке x = p /4 по формуле (22). Итог работы программы приведен на рис.3.

 

Рис.3. Зависимость ошибки аппроксимации второй производной по
формуле (22) в зависимости от шага сетки

 

Листинг_№3

%Программа анализа зависимости ошибки

%аппроксимации второй производной функции

%в зависимости от величины шага сетки

%в фиксированной точке

%очищаем рабочую область

clear all

%выбираем начальное значение шага сетки

h(1)=1.0;

%определяем число делений шага пополам

n=30;

%формируем массив шагов сетки

for i=2:n

h(i)=h(i-1)/2;

end

%выбираем фиксированную точку для

%функции численной оценки второй

%производной функции y=sin(x)

x=pi/4;

%находим численное значение второй производной

%при различных значениях h и сравниваем их с

%эталоном

for i=1:n

d2y=(sin(x+h(i))-2*sin(x)+sin(x-h(i)))/h(i)^2;

er(i)=abs(-sin(x)-d2y);

end

%строим график зависимости ошибки аппроксимации

%в зависимости от величины шага сетки h

plot(h,er);

 

Из графика на рис.3 видно, что ошибка аппроксимации имеет отчетливый минимум, т.е. и очень большое значение шага, и очень малое приводят к заметным, по сравнению с оптимальным значением, ошибкам. Из графика видно, что оптимальное значение шага конечной разности лежит между 10^-4 и 10^-3.

Метод Рунге-Ромберга

Данный метод позволяет получить более высокий порядок точности, не прибегая к более сложным формулам типа (14), (15), в которых учитывается дополнительное по сравнению с минимальным количество узлов аппроксимации. Повышение точности достигается путем комбинирования пары расчетов на разных сетках с разным количеством узлов.

Пример расчета аппроксимации второй производной в (17) позволяет предположить, что и в общем случае, если есть некоторая формула x(x,h) для приближенного вычисления величины g(x), то ошибку аппроксимации можно представить в виде:

                                              (27)

Осуществим еще один расчет по схеме x(x,h) для той же точки с помощью равномерной сетки, но с шагом rh. Для новой сетки ошибка аппроксимации будет иметь вид, аналогичный формуле (27), т.е.

                                     (28)

Имея два расчета (27), (28) для двух сеток можно оценить погрешность R. Она может быть получена после вычитания уравнения (27) из (28), т.е.

.                                   (29)

где считается, что . Формулу (29) принято именовать первой формулой Рунге.

Первая формула Рунге позволяет уточнить погрешность в исходной схеме (27), а именно, подставляя погрешность из (29) в (27), находим

.                                (30)

Формула (30) именуется второй формулой Рунге, и она позволяет получить численный результат с более высоким порядком точности.

Метод Рунге обобщается на случай произвольного количества сеток q. В этом случае, соответствующим образом модифицировав схему расчетов, можно повысить точность аппроксимации до уровня O(h^{p + q-1}). Такая схема расчетов называется схемой Ромберга. Более подробно эта схема изложена в учебнике[1].

Рассмотрим пример, иллюстрирующий работу метода Рунге. Выберем для тестирования схему численного дифференцирования (18), т.е. правую конечную разность, которая, как было установлено выше, имеет первый порядок аппроксимации, т.е. p = 1. Уточнение по Рунге должно повысить точность до второго порядка по шагу сетки. Как и выше, тестируемой функцией будет выступать y =sin(x).

На рис.4 приведена схема позиционированию двух сеток по методу Рунге, когда одна из них вдвое более подробная, чем другая.

 

Рис.4. Схема позиционирования двух сеток в методе Рунге

 

На листинге_№4 приведен код программы, которая уточняет производную синуса по методу Рунге после пары расчетов на исходной сетке и на сетке с удвоенным количеством узлов. Итог работы программы сконцентрирован в графиках рис.5.

 

Листинг_№4

%Программа иллюстрирующая метод

%Рунге по повышению точности численного

%дифференцирования путем комбинирования

%пары расчетов на двух равномерных сетках,

%причем вторая содержит вдвое большее

%количество узлов

%очищаем рабочее пространство

clear all

%определяем шаг исходной сетки

h=0.1;

%формируем исходную сетку

x=0:h:pi;

%определяем число узлов, входящих в исходную

%сетку

n=length(x);

%определяем шаг более подробной сетки

hm=h/2;

%создаем сетку вдвое более подробную

xm=0:hm:pi;

%определяем число узлов, входящих в более

%подробную сетку

m=length(xm);

%рассчитываем производную с помощью правой

%разности на исходной сетке и оцениваем

%соответствующую абсолютную ошибку

for i=1:(n-1)

dy(i)=(sin(x(i+1))-sin(x(i)))/h;

er1(i)=abs(cos(x(i))-dy(i));

end

%рассчитываем производную с помощью правой

%разности на более частой сетке и оцениваем

%соответствующую абсолютную ошибку

for i=1:(m-1)

dym(i)=(sin(xm(i+1))-sin(xm(i)))/hm;

er2(i)=abs(cos(xm(i))-dym(i));

end

%уточняем значение производной с помощью

%метода Рунге

for i=1:(n-1)

dyrunge(i)=dy(i)-2*(dy(i)-dym(2*i-1));

er3(i)=abs(cos(x(i))-dyrunge(i));

end

%строим общий график со всеми тремя кривыми

%ошибок

plot(x([1:(n-1)]),er1([1:(n-1)]),'-o',...

xm([1:2:(m-1)]),er2([1:2:(m-1)]),'-p',...

x([1:(n-1)]),er3([1:(n-1)]),'-h');

 

Сравнение графиков на рис.5 подтверждает теоретические выводы. Процедура Рунге действительно резко повышает точность численной оценки производной в нашем примере.

 

Рис.5. Ошибки двух схем численного дифференцирования: для
исходной схемы и для схемы с удвоенным количеством узлов,
а также ошибка уточняющей процедуры Рунге

 

Рассмотрим еще одну задачу, иллюстрирующую методы численного дифференцирования. Требуется изучить скорость и ускорение динамики народонаселения в Российской Федерации. Данные возьмем из российского статистического ежегодника[2]. Это типичный пример, когда функция задана таблично и требуется найти первую и вторую производную. Эту задачу, естественно, решать в среднем, поскольку значения функции определены с определенными ошибками. Согласно процедуре решения в среднем, неизвестная функция аппроксимируется некоторым полиномом, коэффициенты которого определяются методом наименьших квадратов. Дифференцируя полученный полином один и два раза, находим соответственно скорость и ускорение демографической динамики в РФ.

На листинге_№5 приведен код соответствующей программы.

 

Листинг_№5

%Программа анализа демографической динамики в РФ

%очищаем рабочее пространство

clear all

%вводим данные: год проведения переписи населения

%поделенный для улучшения работы алгоритма на 1000

x=[1.959 1.970 1.979 1.989 1.992 1.993 1.994 1.995 2.002];

%вводим данные: количество народонаселения в млн. чел.

%на соответствующий год переписи населения

y=[117.5 129.9 137.4 147 148.3 148.3 148 147.9 145.2];

%задаем степень аппроксимирующего полинома

nm=3;

%обращаемся к стандартной программе, вычисляющей

%коэффициенты полиномы, аппроксимирующего данные

%в среднем

p=polyfit(x,y,nm);

%определяем значения аппроксимирующего полинома

%в отчетные моменты времени

phi=polyval(p,x);

%рисуем исходные данные совместно с аппроксимирующим

%полиномом

plot(1000*x,y,'*',1000*x,phi);

%определяем коэффициенты полинома, описывающего

%скорость демографической динамики

for i=1:nm

dp1(i)=(nm-i+1)*p(i);

end

 

Рис.6. График демографической динамики в РФ

 

%определяем значения аппроксимирующего полинома

%скорости демографической динамики в отчетные

%моменты времени

dphi1=polyval(dp1,x);

%рисуем график скорости демографической динамики

%plot(1000*x,dphi1/1000);

%определяем значения аппроксимирующего полинома

%ускорения демографической динамики в отчетные

%моменты времени

for i=1:(nm-1)

dp2(i)=(nm-i+1)*(nm-i)*p(i);

end

%определяем значения аппроксимирующего полинома

%ускорения демографической динамики в отчетные

%моменты времени

dphi2=polyval(dp2,x);

%рисуем график ускорения демографической динамики

%plot(1000*x,dphi2/1e6);

grid on

 

На рис.6 приведен график демографической динамики в РФ. Маркерами отмечены табличные значения, линия соответствует кубической параболе наилучшим образом аппроксимирующей наши данные в смысле метода наименьших квадратов. Из графика на рис.6 видно, что народонаселение в РФ, начиная с середины 90-х годов, неуклонно сокращается.

На рис.7,а приведен график скорости демографической динамики, а на рис.7,б — график ускорения.

 

Рис.7,а. Демографическая скорость в РФ	Рис.7,б. Демографическое ускорение в РФ

 

Согласно графику на рис.7,а, скорость демографической динамики сменила знак и стала отрицательной в середине 90-х годов, что согласуется с графиком демографической динамики на рис.6. Ускорение демографической динамики сменило знак и стало отрицательным в начале 70-х годов, что также согласуется с графиком демографической динамики, где на начало 70-х годов приходится пик скорости прироста населения.

[1] Калиткин Н.Н. Численные методы. — М.: Гл. ред. физ.-мат лит., 1978. 512с.

[2] Российский статистический ежегодник. Стат. Сб./ Госкомстат России. М., 2001.

---

Дата добавления: 2022-01-22; просмотров: 36; Мы поможем в написании вашей работы!

Поделиться с друзьями:

---

---

Мы поможем в написании ваших работ!