Применение интерполирования при дифференцировании

Материал из MachineLearning.

Перейти к: навигация, поиск

Содержание

Введение

Постановка математической задачи

Численное дифференцирование применяется, если функцию y(x) трудно или невозможно продифференцировать аналитически - например, если она задана таблицей. Оно нужно также при решении дифференциальных уравнений при помощи разностных методов.

Изложение метода

При численном дифференцировании функцию y(x) аппроксимируют легко вычисляемой функцией \varphi(x) и приближенно полагают y'(x)=\varphi'(x). При этом можно использовать различные способы аппроксимации.

Интерполирование полиномами Ньютона

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

Пусть задана сетка x_0<x_1<\dots<x_N. y(x) - исследуемая функция. Введем обозначение \xi_i=x-x_i и введем понятие разделенные разности y(x_0,x_1,\dots,x_k), определяемые следующим образом:

y(x_i,x_j)=\frac{y(x_i) - y(x_j)}{x_i - x_j},
y(x_i,x_j,x_k) = \frac{y(x_i,x_j)-y(x_j,x_k)}{x_i-x_k} и т.д.

Можно доказать, что

y(x_0,x_1,\dots,x_k) = \sum_{p=0}^{k}y_p \prod_{i=0, i\neq p}^k {(x_p-x_i)}^{-1}.

Запишем интерполяционный многочлен Ньютона и продифференцируем его почленно:

\varphi(x)=y(x_0)+\xi_0 y(x_0,x_1)+\xi_0\xi_1 y(x_0,x_1,x_2) + \xi_0\xi_1\xi_2  y(x_0,x_1,x_2,x_3)+\dots
\varphi'(x)=y(x_0,x_1)+(\xi_0+\xi_1) y(x_0,x_1,x_2) + (\xi_0\xi_1+\xi_0\xi_2 +\xi_1\xi_2) y(x_0,x_1,x_2,x_3)+\dots
\varphi''(x)=2 y(x_0,x_1,x_2) + 2(\xi_0+\xi_1 +\xi_2) y(x_0,x_1,x_2,x_3)+\dots

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

( 1 )
 \varphi^{(k)}(x)=k!\left[ y(x_0,x_1,\dots,x_k) + \left( \sum_{i=0}^k \xi_i \right) y(x_0,x_1,\dots,x_{k+1}) + \left( \sum_{i>j\geq 0}^{i=k+1}\xi_i\xi_j\right)y(x_0,x_1,\dots,x_{k+2})   + \left( \sum_{i>j>l\geq 0}^{i=k+2}\xi_i\xi_j\xi_l\right)y(x_0,x_1,\dots,x_{k+3}) +\dots\right]

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

( 2 )
y'(x)\approx y(x_0,x_1) = \frac{y(x_0)-y(x_1)}{x_0-x_1},
\frac{1}{2}y''(x)\approx y(x_0,x_1,x_2) = \frac{1}{x_0-x_2}\left( \frac{y_0-y_1}{x_0-x_1}- \frac{y_1-y_2}{x_1-x_2}\right),
\frac{1}{k!} y^{(k)}(x) \approx y(x_0,x_1,\dots,x_k) = \sum_{p=0}^{k}y_p \prod_{i=0, i\neq p}^k {(x_p-x_i)}^{-1}

Исследование точности полученных выражений при численных расчётах удобно делать при помощи апостериорной оценки, по скорости убывания членов ряда (1). Если шаг сетки достаточно мал, то погрешность близка к первому отброшенному члену. Пусть мы используем узлы x_i, i=1\dots n. Тогда первый отброшенный член содержит разделенную разность y(x_0,x_1,\dots,x_{n+1}), которая согласно (2) примерно равна y^{(n+1)}(x)/(n+1)!. Перед ней стоит сумма произведений различных множителей \xi_i ; каждое произведение содержит n+1-k множителей, а вся сумма состоит из C_{n+1}^k слагаемых. Отсюда следует оценка погрешности формулы (3) с  n+1 узлами:

(3)
R_n^{(k)}\leq\frac{M_{n+1}}{(n+1-k)!}\max_i {|\xi_i|}^{n+1-k},\; M_{n+1}=\max{|y^{(n+1)}|}

В частности, если сетка равномерная, то M_{n+1}=\max{|\xi_i|}\leq nh, откуда

(4)
R_n^{(k)}<M_{n+1} {\left( \frac{en}{n+1-k}h \right)}^{n+1-k}=O(h^{n+1-k}) .

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

Таким образом, порядок точности формулы (1) по отношению к шагу сетки равен числу оставленных в ней членов, или, что то же самое, он равен числу узлов интерполяции минус порядок производной. Поэтому минимальное число узлов, необходимое для вычисления k-й производной, равно k+1; оно приводит к формулам (2) и обеспечивает первый порядок точности. Эти выводы соответствуют общему принципу: при почленном дифференцировании ряда скорость его сходимости уменьшается.

Интерполирование полиномами Лагранжа

Рассмотрим неравномерную сетку \omega_h=\{a=x_0<x_1<x_2<\dots<x_N=b\} и обозначим за h_i=x_i-x_{i-1}, i=1,2,\dots,N шаги этой сетки. В качества примера получим формулы численного дифференцирования, основанные на использовании многочлена Лагранжа L_{2,i}(x), построенного для функции u(x) по трем точкам x_{i-1},x_i,x_{i+1}. Многочлен L_{2,i}(x) имеет вид

L_{2,i}(x)=\frac {(x-x_i)(x-x_{i+1})}{h_i(h_i+h_{i+1})}u_{i-1}-\frac {(x-x_{i-1})(x-x_{i+1})}{h_ih_{i+1}}u_i+\frac {(x-x_{i-1})(x-x_i)}{h_{i+1}(h_i+h_{i+1})}u_{i+1}

Отсюда получим L_{2,i}'(x)=\frac{(2x-x_i-x_{i+1})}{h_i(h_i+h_{i+1})}u_{i-1}-\frac {(2x-x_{i-1}-x_{i+1})}{h_ih_{i+1}}u_i+\frac {(2x-x_{i-1}-x_i)}{h_{i+1}(h_i+h_{i+1})}u_{i+1}

Это выражение можно принять за приближенное значение u'(x) в любой точке x[x_{i-1},x_{i+1}]. Его удобнее записать в виде  L_{2,i}'(x)=\frac {1}{\bar{h_i}}[(x-x_{i-\frac{1}{2}}) \frac{u_{i+1}-u_i}{h_{i+1}} + (x_{i+\frac{1}{2}}-x) \frac{u_i-u_{i-1}}{h_i}] , где \bar{h_i}=0,5(h_i+h_{i+1}), x_{i-\frac{1}{2}}=x_i-0,5h_i.

В частности, при x=x_i получим L_{2,i}'(x_i)=\frac{1}{2}(\frac{h_i}{\bar{h_i}}\frac{u_{i+1}-u_i}{h_{i+1}}+\frac{h_{i+1}}{\bar{h_i}}\frac{u_i-u_{i-1}}{h_i}), И если сетка равномерна, h_{i+1}=h_i=h, то приходим к центральной разностной производной, L_{2,i}'(x_i)=u_{\dot{x},i}. При использовании интерполяционного многочлена первой степени точно таким образом можно получить односторонние разностные производные u_{\bar{x},i} и u_{x,i}. Далее вычисляя вторую производную многочлена L_{2,i}(x), получим приближенное выражение для u''(x) при x[x_{i-1},x_{i+1}]:

u''(x)L_{2,i}''(x)=\frac{1}{\bar{h_i}}(\frac{u_{i+1}-u_i}{h_{i+1}}- \frac{u_i-u_{i-1}}{h_i})

На равномерной сетке это выражение совпадает со второй разностной производной u_{\bar{x}x,i}. Ясно, что для приближенного вычисления дальнейших производных уже недостаточно многочлена L_{2,i}(x), надо привлекать многочлены более высокого порядка и тем самым увеличивать число узлов, участвующих в аппроксимации.

Порядок погрешности аппроксимации зависит как от порядка интерполяционного многочлена, так и от расположения узлов интерполирования. Получим выражение для погрешности аппроксимации, возникающей при замене u'(x) выражением L_{2,i}'(x). Будем считать, что x[x_{i-1},x_{i+1}] и что величины h_i, h_{i+1} имеют один и тот же порядок малости при измельчении сетки. По формуле Тейлора в предположении ограниченности u^{(4)}(x) получим u_{i+k}=u(x)+(x_{i+k}-x)u'(x)+\frac{{(x_{i+k}-x)}^2}{2}u''(x) +\frac{{(x_{i+k}-x)}^3}{3}u'''(x) +O(h^4),

где k=01,h=max\{h_i,h_{i+1}\}

Отсюда приходим к следующим разложениям разностных отношений

\frac{u_i-u_{i-1}}{h_i}=u'(x)-(x-x_{i-\frac{1}{2}})u''(x)+(\frac{{(x-x_{i-\frac{1}{2}})}^2}{2}+\frac{h_i^2}{24})u'''(x)+O(h^3)

\frac{u_{i+1}-u_i}{h_{i+1}}=u'(x)-(x_{i+\frac{1}{2}}-x)u''(x)+(\frac{{(x_{i+\frac{1}{2}}-x)}^2}{2}+\frac{h_{i+1}^2}{24})u'''(x)+O(h^3)

Подставляя полученные формулы в выражение для разностной производной и приводя подобные слагаемые получим

L_{2,i}'(x)=u'(x)-[\frac{{(x-x_i)}^2}{2}-\frac{(h_{i+1}-h_i)(x-x_i)}{3}-\frac{h_ih_{i+1}}{6}]u'''(x)+O(h^3), x[x_{i-1},x_{i+1}].

Отсюда видно,что разностное выражение аппроксимирует u'(x) со вторым порядком.

Если подставить полученные ранее разностные отношения в выражение для второй производной многочлена L_{2,i}(x), то имеем

L_{2,i}''(x)=u''(x)+(x_i-x + \frac{h_{i+1}-h_i}{3})u'''(x)+O(h^2)

Из этого выражения видно, что даже на равномерной сетке,т.е. когда h_i=h_{i+1}, второй порядок аппроксимации имеет место лишь в точке x=x_i, а относительно других точек (например,x=x_{i+1}) выполняется аппроксимация только первого порядка. Таким образом, получим аппроксимацию лишь первого порядка.

Интеполирование кубическими сплайнами

Для того, чтобы избежать больших погрешностей в процессе приближения, весь отрезок [a,b] разбивают на частичные отрезки и на каждом из частичных отрезков приближенно заменяют функцию f(x) многочленом невысокой степени (так называемая кусочно-полиномиальная интерполяция). Одним из способов интерполирования на всем отрезке является интерполирование с помощью сплайн-функций.

Сплайн-функцией или сплайном называют кусочно-полиномиальную функцию, определенную на отрезке [a,b] и имеющую на этом отрезке некоторое число непрерывных производных. Преимущество сплайнов перед обычной интерполяцией является, во-первых их сходимость, во-вторых, устойчивость процесса вычислений.

Построение кубического сплайна.

Пусть на [a,b] задана непрерывная функция f(x). Введем сетку a=x_0<x_1<x_2<\dots<x_N=b и обозначим f_i=f(x_i), i=0,1,\dots,N. Сплайном соответствующим данной функции f(x) и данным узлам \{x_i\}_{i=0}^N называется функция s(x), удовлетворяющая следующим условиям:

а) на каждом сегменте [x_i-1,x_i], i=1,2,\dots,N функция s(x) является многочленом третьей степени;

б) функция s(x), а также её первая и вторая производная производные непрерывны на [a,b];

в) s(x_i)=f(x_i),i=1,2,\dots,N;

На каждом из отрезков [x_i-1,x_i], i=1,2,\dots,N будем искать функцию s(x)=s_i(x) в виде многочлена третьей степени

s_i(x)=a_i+b_i(x-x_i)+\frac{c_i}{2}(x-x_i)^2 +\frac{d_i}{6}(x-x_i)^3,

где x_i-1<=x<=x_i, i=1,2,\dots,N, где a_i,b_i,c_i,d_i - коэффициенты, подлежащие определению. Доказано, что существует единственный кубический сплайн, определяемый условиями а)-в) и граничными условиями s''(a)=s''(b)=0.

Для их нахождения используются следующие формулы

1) a_i=f(x_i),i=1,2,\dots,N

2) Для определения коэффициентов c_i получаем систему уравнений

h_ic_{i-1}+2(h_i+h_{i+1})c_i+h_{i+1}c_{i+1}=6\left(\frac{f_{i+1}-f_i}{h_{i+1}}-\frac{f_i-f_{i-1}}{h_i}\right), <tex>i=1,2,\dots,N-1 ,c_0=c_N=0

(система решается методом прогонки) По найденным коэффициентам c_i коэффициенты b_i, d_i определяются с помощью явных формул

3) d_i=\frac{c_i-c_{i-1}}{h_i},<tex>i=1,2,\dots,N

4) b_i=\frac{h_i}{2}c_i-\frac{h_i^2}{6}d_i+\frac{f_i-f_{i-1}}{h_i},
<tex>i=1,2,\dots,N

Найдем производные введенного кубического сплайна, имеем s_i'(x)=b_i+c_i(x-x_i)+\frac{d_i}{2}{(x-x_i)}^2

s_i''(x)=c_i+d_i(x-x_i)

s_i'''(x)=d_i

Рассмотрим оценку погрешности метода, которая зависит от выбора сеток и от гладкости f(x). Для простоты изложения допустим, что сетка равномерная, т.е.

\omega_h=\{x_i=a+ih, i=0,1,\dots,N\} с шагом b=\frac{b-a}{N}

От функции f(x) будем требовать существования непрерывной на [a,b] четвертой производной, f(x)C^{(4)}[a,b]. Кроме того, предположим, что выполнены граничные условия f''(a)=f''(b)=0 и такие же условия для сплайнов. Обозначим,

||g(x)||_{C[a,b]}=\max_{[a,b]}|g(x)|, M_4=||f^4(x)||_{C[a,b]}

Пусть s_h(x) - кубический сплайн, построенный для функции f(x) на сетке \omega_h. В следующей теореме приведены оценки погрешности интерполяции для функции f(x) и её производных f'(x), f''(x)

Теорема

Для f(x)C^{(4)}[a,b] справедливы оценки

||f(x)-s_h(x)||_{C[a,b]}<=M_4h^4

||f'(x)-s_h'(x)||_{C[a,b]}<=M_4h^3

||f''(x)-s_h''(x)||_{C[a,b]}<=M_4h^2

Из этих оценок следует, что при h \to 0 (т.е. при N \to \infty) последовательности s_h^{(i)}(x), i=0,1,2 сходятся соответственно к функциям f^{(i)}(x) i=0,1,2.

Обычно дифференцирование кубического сплайна позволят определить первую и вторую производную интерполяционного многочлена с хорошей точностью. Если надо вычислить более высокие производные, то целесообразно строить сплайны высоких порядков. Из-за большей трудоемкости этот способ редко используется. Способ дифференцирования с помощью сплайновой интерполяцией теоретически мало исследован.

Тригонометрическая интерполяция

Не всякую функцию целесообразно приближать алгебраическими многочленами. Рассмотрим тригонометрическую интерполяцию. Если f(x) - периодическая функция с периодом l, то естественно строить приближения с помощью функций

\varphi_k(x)=a_k\cos(\frac{\pi kx}{l})+b_k\sin(\frac{\pi kx}{l}), k=0,1,\dots,n

Таким образом, тригонометрическая интерполяция состоит в замене f(x) тригонометрическим многочленом

T_n(x)=\sum_{k=0}^n \varphi_k(x)= a_0 + \sum_{k=1}^n\left(a_k\cos(\frac{\pi kx}{l})+b_k\sin(\frac{\p ikx}{l})\right),

коэффициенты которого отыскиваются из системы уравнений

T_n(x_j)=f(x_j), j=1,2,\dots,2n+1,

где x_0<x_1<\dots<x_{2n+1},x_{2n+1}-x_0=l.

Рекомендации программисту

При реализации достаточно сложно и трудоемко использовать методы сплайнов и метод тригонометрической интерполяции. Поэтому был рассмотрен случай интерполирования полиномом Лагранжа по трем точкам. При вычислении полиномов ошибка была ничтожна. При попытке в качестве функции взять синус,ошибка достигала очень большх цифр и был порядка 1.0, что достаточно велико для этой функции. Поэтому следует внимательно выбирать шаг, расстояние между точками. Программа требует ввода координат трех точек для полинома x^5 в порядке возрастания,а также точки промежуточной, в которой будет находиться значение производной. На экран выводится значение, полученное при помощи полинома Лагранжа,истинное значение, полученное аналитически и ошибка вычисления. Приведем примеры работы программы

Координаты первой точки x_1 Координаты второй точки x_2 Координаты третьей точки x_3 Координаты четвертой точки x_4 Ошибка вычисления
0.1 0.2 0.3 0.15 0.0056875
1.0 2.0 3.0 1.5 5.6875
0.01 0.02 0.03 0.015 -0.00250664
0.15 0.24 0.3 0.22 0.0042998

Список литературы

  • А.А.Самарский, А.В.Гулин.  Численные методы. Москва «Наука», 1989.
  • Н.Н.Калиткин.  Численные методы. Москва «Наука», 1978.

См. также

Личные инструменты