Участник:Айнагуль Джумабекова/Песочница
Материал из MachineLearning.
(9 промежуточных версий не показаны.) | |||
Строка 1: | Строка 1: | ||
- | <tex> | + | == Введение == |
+ | === Постановка математической задачи === | ||
+ | Численное дифференцирование применяется, если функцию <tex>y(x)</tex> трудно или невозможно продифференцировать аналитически - например, если она задана таблицей. Оно нужно также при решении дифференциальных уравнений при помощи разностных методов. | ||
- | + | == Изложение метода == | |
- | <tex> | + | При численном дифференцировании функцию <tex>y(x)</tex> аппроксимируют легко вычисляемой функцией <tex>\varphi(x)</tex> и приближенно полагают |
+ | <tex>y'(x)=\varphi'(x)</tex>. При этом можно использовать различные способы аппроксимации. | ||
- | + | ===Интерполирование полиномами Ньютона=== | |
+ | Рассмотрим случай аппроксимации интерполяционным многочленом Ньютона. | ||
- | <tex>L_{2,i}'(x_i)=u_{\dot{x},i}</tex> | + | Пусть задана сетка <tex>x_0<x_1<\dots<x_N. y(x)</tex> - исследуемая функция. Введем обозначение <tex>\xi_i=x-x_i</tex> и введем понятие разделенные разности <tex>y(x_0,x_1,\dots,x_k)</tex>, определяемые следующим образом: |
+ | |||
+ | ::<tex>y(x_i,x_j)=\frac{y(x_i) - y(x_j)}{x_i - x_j}</tex>, | ||
+ | |||
+ | ::<tex>y(x_i,x_j,x_k) = \frac{y(x_i,x_j)-y(x_j,x_k)}{x_i-x_k}</tex> и т.д. | ||
+ | |||
+ | Можно доказать, что | ||
+ | |||
+ | ::<tex>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}</tex>. | ||
+ | |||
+ | Запишем интерполяционный многочлен Ньютона и продифференцируем его почленно: | ||
+ | |||
+ | ::<tex>\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</tex> | ||
+ | |||
+ | ::<tex>\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</tex> | ||
+ | |||
+ | ::<tex>\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</tex> | ||
+ | |||
+ | Общая формула примет следующий вид: | ||
+ | {{ eqno | 1 }} | ||
+ | ::<tex> \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] </tex> | ||
+ | |||
+ | Обрывая ряд на некотором числе членов, получим приближенное выражение для соответствующей производной. Наиболее простые выражения получим, оставляя в формуле {{eqref|1}} только первый член: | ||
+ | {{ eqno | 2 }} | ||
+ | ::<tex>y'(x)\approx y(x_0,x_1) = \frac{y(x_0)-y(x_1)}{x_0-x_1}</tex>, | ||
+ | |||
+ | ::<tex>\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)</tex>, | ||
+ | |||
+ | ::<tex>\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}</tex> | ||
+ | |||
+ | Исследование точности полученных выражений при численных расчётах удобно делать при помощи апостериорной оценки, по скорости убывания членов ряда {{eqref|1}}. Если шаг сетки достаточно мал, то погрешность близка к первому отброшенному члену. Пусть мы используем узлы <tex>x_i, i=1\dots n</tex>. Тогда первый отброшенный член содержит разделенную разность <tex>y(x_0,x_1,\dots,x_{n+1})</tex>, которая согласно {{eqref|2}} примерно равна <tex>y^{(n+1)}(x)/(n+1)!</tex>. Перед ней стоит сумма произведений различных множителей <tex>\xi_i </tex>; каждое произведение содержит <tex>n+1-k</tex> множителей, а вся сумма состоит из <tex>C_{n+1}^k</tex> слагаемых. Отсюда следует оценка погрешности формулы {{eqref|3}} с <tex> n+1 </tex> узлами: | ||
+ | |||
+ | {{eqno|3}} | ||
+ | ::<tex>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)}|}</tex> | ||
+ | |||
+ | В частности, если сетка равномерная, то <tex>M_{n+1}=\max{|\xi_i|}\leq nh</tex>, откуда | ||
+ | {{eqno|4}} | ||
+ | ::<tex>R_n^{(k)}<M_{n+1} {\left( \frac{en}{n+1-k}h \right)}^{n+1-k}=O(h^{n+1-k}) </tex>. | ||
+ | |||
+ | Стоит заметить, что строгое априорное исследование погрешности формулы {{eqref|3}}, аналогичное выводу остаточного члена многочлена Ньютона в форме Коши, для произвольного расположения узлов приводит к той же оценке {{eqref|3}}. | ||
+ | |||
+ | Таким образом, порядок точности формулы {{eqref|1}} по отношению к шагу сетки равен числу оставленных в ней членов, или, что то же самое, он равен числу узлов интерполяции минус порядок производной. Поэтому минимальное число узлов, необходимое для вычисления <tex>k</tex>-й производной, равно <tex>k+1</tex>; оно приводит к формулам {{eqref|2}} и обеспечивает первый порядок точности. Эти выводы соответствуют общему принципу: при почленном дифференцировании ряда скорость его сходимости уменьшается. | ||
+ | |||
+ | ===Интерполирование полиномами Лагранжа=== | ||
+ | |||
+ | Рассмотрим неравномерную сетку <tex>\omega_h=\{a=x_0<x_1<x_2<\dots<x_N=b\}</tex> | ||
+ | и обозначим за <tex>h_i=x_i-x_{i-1}</tex>, <tex>i=1,2,\dots,N</tex> шаги этой сетки. В качества примера получим формулы численного дифференцирования, основанные на использовании многочлена Лагранжа <tex>L_{2,i}(x)</tex>, построенного для функции <tex>u(x)</tex> по трем точкам <tex>x_{i-1},x_i,x_{i+1}</tex>. | ||
+ | Многочлен <tex>L_{2,i}(x)</tex> имеет вид | ||
+ | |||
+ | <tex>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}</tex> | ||
+ | |||
+ | Отсюда получим | ||
+ | <tex>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}</tex> | ||
+ | |||
+ | Это выражение можно принять за приближенное значение <tex>u'(x)</tex> в любой точке <tex>x</tex>∈ <tex>[x_{i-1},x_{i+1}]</tex>. | ||
+ | Его удобнее записать в виде | ||
+ | <tex> 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}]</tex> , где | ||
+ | <tex>\bar{h_i}=0,5(h_i+h_{i+1})</tex>, <tex>x_{i-\frac{1}{2}}=x_i-0,5h_i</tex>. | ||
+ | |||
+ | В частности, при <tex>x=x_i</tex> получим | ||
+ | <tex>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})</tex>, | ||
+ | И если сетка равномерна, <tex>h_{i+1}=h_i=h</tex>, то приходим к центральной разностной производной, <tex>L_{2,i}'(x_i)=u_{\dot{x},i}</tex>. | ||
+ | При использовании интерполяционного многочлена первой степени точно таким образом можно получить односторонние разностные производные <tex>u_{\bar{x},i}</tex> и <tex>u_{x,i}</tex>. | ||
+ | Далее вычисляя вторую производную многочлена <tex>L_{2,i}(x)</tex>, получим приближенное выражение для <tex>u''(x)</tex> при <tex>x</tex>∈<tex>[x_{i-1},x_{i+1}]</tex>: | ||
<tex>u''(x)</tex>≈<tex>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})</tex> | <tex>u''(x)</tex>≈<tex>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})</tex> | ||
+ | |||
+ | На равномерной сетке это выражение совпадает со второй разностной производной <tex>u_{\bar{x}x,i}</tex>. Ясно, что для приближенного вычисления дальнейших производных уже недостаточно многочлена <tex>L_{2,i}(x)</tex>, надо привлекать многочлены более высокого порядка и тем самым увеличивать число узлов, участвующих в аппроксимации. | ||
+ | |||
+ | Порядок погрешности аппроксимации зависит как от порядка интерполяционного многочлена, так и от расположения узлов интерполирования. Получим выражение для погрешности аппроксимации, возникающей при замене <tex>u'(x)</tex> выражением <tex>L_{2,i}'(x)</tex>. Будем считать, что <tex>x</tex>∈ <tex>[x_{i-1},x_{i+1}]</tex> и что величины <tex>h_i, h_{i+1}</tex> имеют один и тот же порядок малости при измельчении сетки. По формуле Тейлора в предположении ограниченности <tex>u^{(4)}(x)</tex> получим | ||
+ | <tex>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)</tex>, | ||
+ | |||
+ | где <tex>k=0</tex>,±<tex>1,h=max\{h_i,h_{i+1}\}</tex> | ||
+ | |||
+ | Отсюда приходим к следующим разложениям разностных отношений | ||
+ | |||
+ | <tex>\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)</tex> | ||
+ | |||
+ | <tex>\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)</tex> | ||
+ | |||
+ | Подставляя полученные формулы в выражение для разностной производной и приводя подобные слагаемые получим | ||
+ | |||
+ | <tex>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)</tex>, <tex>x</tex>∈ <tex>[x_{i-1},x_{i+1}]</tex>. | ||
+ | |||
+ | Отсюда видно,что разностное выражение аппроксимирует <tex>u'(x)</tex> со вторым порядком. | ||
+ | |||
+ | Если подставить полученные ранее разностные отношения в выражение для второй производной многочлена <tex>L_{2,i}(x)</tex>, то имеем | ||
+ | |||
+ | <tex>L_{2,i}''(x)=u''(x)+(x_i-x + \frac{h_{i+1}-h_i}{3})u'''(x)+O(h^2)</tex> | ||
+ | |||
+ | Из этого выражения видно, что даже на равномерной сетке,т.е. когда <tex>h_i=h_{i+1}</tex>, второй порядок аппроксимации имеет место лишь в точке <tex>x=x_i</tex>, а относительно других точек (например,<tex>x=x_{i+1}</tex>) выполняется аппроксимация только первого порядка. | ||
+ | Таким образом, получим аппроксимацию лишь первого порядка. | ||
+ | |||
+ | ===Интеполирование кубическими сплайнами=== | ||
+ | Для того, чтобы избежать больших погрешностей в процессе приближения, весь отрезок [a,b] разбивают на частичные отрезки и на каждом из частичных отрезков приближенно заменяют функцию <tex>f(x)</tex> многочленом невысокой степени (так называемая кусочно-полиномиальная интерполяция). | ||
+ | Одним из способов интерполирования на всем отрезке является интерполирование с помощью сплайн-функций. | ||
+ | |||
+ | Сплайн-функцией или сплайном называют кусочно-полиномиальную функцию, определенную на отрезке [a,b] и имеющую на этом отрезке некоторое число непрерывных производных. | ||
+ | Преимущество сплайнов перед обычной интерполяцией является, во-первых их сходимость, во-вторых, устойчивость процесса вычислений. | ||
+ | |||
+ | Построение кубического сплайна. | ||
+ | |||
+ | Пусть на [a,b] задана непрерывная функция <tex>f(x)</tex>. Введем сетку | ||
+ | <tex>a=x_0<x_1<x_2<\dots<x_N=b</tex> | ||
+ | и обозначим <tex>f_i=f(x_i)</tex>, <tex>i=0,1,\dots,N</tex>. | ||
+ | Сплайном соответствующим данной функции <tex>f(x)</tex> и данным узлам <tex>\{x_i\}_{i=0}^N</tex> называется функция <tex>s(x)</tex>, удовлетворяющая следующим условиям: | ||
+ | |||
+ | а) на каждом сегменте <tex>[x_i-1,x_i]</tex>, <tex>i=1,2,\dots,N</tex> функция <tex>s(x)</tex> является многочленом третьей степени; | ||
+ | |||
+ | б) функция <tex>s(x)</tex>, а также её первая и вторая производная производные непрерывны на [a,b]; | ||
+ | |||
+ | в) <tex>s(x_i)=f(x_i)</tex>,<tex>i=1,2,\dots,N</tex>; | ||
+ | |||
+ | На каждом из отрезков <tex>[x_i-1,x_i]</tex>, <tex>i=1,2,\dots,N</tex> будем искать функцию <tex>s(x)=s_i(x)</tex> в виде многочлена третьей степени | ||
+ | |||
+ | <tex>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</tex>, | ||
+ | |||
+ | где <tex>x_i-1<=x<=x_i</tex>, <tex>i=1,2,\dots,N</tex>, где <tex>a_i,b_i,c_i,d_i</tex> - коэффициенты, подлежащие определению. | ||
+ | Доказано, что существует единственный кубический сплайн, определяемый условиями а)-в) и граничными условиями <tex>s''(a)=s''(b)=0</tex>. | ||
+ | |||
+ | Для их нахождения используются следующие формулы | ||
+ | |||
+ | 1) <tex>a_i=f(x_i)</tex>,<tex>i=1,2,\dots,N</tex> | ||
+ | |||
+ | 2) Для определения коэффициентов <tex>c_i</tex> получаем систему уравнений | ||
+ | |||
+ | <tex>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</tex> ,<tex>c_0=c_N=0</tex> | ||
+ | |||
+ | (система решается методом прогонки) | ||
+ | По найденным коэффициентам <tex>c_i</tex> коэффициенты <tex>b_i</tex>, <tex>d_i</tex> определяются с помощью явных формул | ||
+ | |||
+ | 3) <tex>d_i=\frac{c_i-c_{i-1}}{h_i},<tex>i=1,2,\dots,N</tex> | ||
+ | |||
+ | 4) <tex>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</tex> | ||
+ | |||
+ | Найдем производные введенного кубического сплайна, имеем | ||
+ | <tex>s_i'(x)=b_i+c_i(x-x_i)+\frac{d_i}{2}{(x-x_i)}^2</tex> | ||
+ | |||
+ | <tex>s_i''(x)=c_i+d_i(x-x_i)</tex> | ||
+ | |||
+ | <tex>s_i'''(x)=d_i</tex> | ||
+ | |||
+ | Рассмотрим оценку погрешности метода, которая зависит от выбора сеток и от гладкости <tex>f(x)</tex>. Для простоты изложения допустим, что сетка равномерная, т.е. | ||
+ | |||
+ | <tex>\omega_h=\{x_i=a+ih, i=0,1,\dots,N\}</tex> с шагом <tex>b=\frac{b-a}{N}</tex> | ||
+ | |||
+ | От функции <tex>f(x)</tex> будем требовать существования непрерывной на [a,b] четвертой производной, <tex>f(x)</tex>∈ <tex>C^{(4)}[a,b]</tex>. Кроме того, предположим, что выполнены граничные условия <tex>f''(a)=f''(b)=0</tex> и такие же условия для сплайнов. Обозначим, | ||
+ | |||
+ | <tex>||g(x)||_{C[a,b]}=\max_{[a,b]}|g(x)|</tex>, <tex>M_4=||f^4(x)||_{C[a,b]}</tex> | ||
+ | |||
+ | Пусть <tex>s_h(x)</tex> - кубический сплайн, построенный для функции <tex>f(x)</tex> на сетке <tex>\omega_h</tex>. В следующей теореме приведены оценки погрешности интерполяции для функции <tex>f(x)</tex> и её производных <tex>f'(x)</tex>, <tex>f''(x)</tex> | ||
+ | |||
+ | '''Теорема''' | ||
+ | |||
+ | Для <tex>f(x)</tex>∈ <tex>C^{(4)}[a,b]</tex> справедливы оценки | ||
+ | |||
+ | <tex>||f(x)-s_h(x)||_{C[a,b]}<=M_4h^4</tex> | ||
+ | |||
+ | <tex>||f'(x)-s_h'(x)||_{C[a,b]}<=M_4h^3</tex> | ||
+ | |||
+ | <tex>||f''(x)-s_h''(x)||_{C[a,b]}<=M_4h^2</tex> | ||
+ | |||
+ | Из этих оценок следует, что при <tex>h \to 0</tex> (т.е. при <tex>N \to \infty</tex>) последовательности <tex>s_h^{(i)}(x)</tex>, <tex>i=0,1,2</tex> сходятся соответственно к функциям <tex>f^{(i)}(x)</tex> <tex>i=0,1,2</tex>. | ||
+ | |||
+ | Обычно дифференцирование кубического сплайна позволят определить первую и вторую производную интерполяционного многочлена с хорошей точностью. Если надо вычислить более высокие производные, то целесообразно строить сплайны высоких порядков. Из-за большей трудоемкости этот способ редко используется. Способ дифференцирования с помощью сплайновой интерполяцией теоретически мало исследован. | ||
+ | |||
+ | ===Тригонометрическая интерполяция=== | ||
+ | Не всякую функцию целесообразно приближать алгебраическими многочленами. | ||
+ | Рассмотрим тригонометрическую интерполяцию. | ||
+ | Если <tex>f(x)</tex> - периодическая функция с периодом l, то естественно строить приближения с помощью функций | ||
+ | |||
+ | <tex>\varphi_k(x)=a_k\cos(\frac{\pi kx}{l})+b_k\sin(\frac{\pi kx}{l}), k=0,1,\dots,n</tex> | ||
+ | |||
+ | Таким образом, тригонометрическая интерполяция состоит в замене <tex>f(x)</tex> тригонометрическим многочленом | ||
+ | |||
+ | ::<tex>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)</tex>, | ||
+ | |||
+ | коэффициенты которого отыскиваются из системы уравнений | ||
+ | |||
+ | <tex>T_n(x_j)=f(x_j), j=1,2,\dots,2n+1</tex>, | ||
+ | |||
+ | где <tex>x_0<x_1<\dots<x_{2n+1},x_{2n+1}-x_0=l</tex>. | ||
+ | |||
+ | ==Рекомендации программисту== | ||
+ | При реализации достаточно сложно и трудоемко использовать методы сплайнов и метод тригонометрической интерполяции. Поэтому был рассмотрен случай интерполирования полиномом Лагранжа по трем точкам. | ||
+ | При вычислении полиномов ошибка была ничтожна. | ||
+ | При попытке в качестве функции взять синус,ошибка достигала очень большх цифр и был порядка 1.0, что достаточно велико для этой функции. | ||
+ | Поэтому следует внимательно выбирать шаг, расстояние между точками. | ||
+ | Программа требует ввода координат трех точек для полинома <tex>x^5</tex> в порядке возрастания,а также точки промежуточной, в которой будет находиться значение производной. | ||
+ | На экран выводится значение, полученное при помощи полинома Лагранжа,истинное значение, полученное аналитически и ошибка вычисления. | ||
+ | |||
+ | == Список литературы == | ||
+ | * ''А.А.Самарский, А.В.Гулин.'' Численные методы. Москва «Наука», 1989. | ||
+ | * ''Н.Н.Калиткин.'' Численные методы. Москва «Наука», 1978. | ||
+ | |||
+ | == См. также == | ||
+ | *[[Практикум ММП ВМК, 4й курс, осень 2008]] | ||
+ | *[[Тригонометрическая интерполяция]] | ||
+ | *[[Интерполяция полиномами Лагранжа и Ньютона]] |
Текущая версия
Содержание |
Введение
Постановка математической задачи
Численное дифференцирование применяется, если функцию трудно или невозможно продифференцировать аналитически - например, если она задана таблицей. Оно нужно также при решении дифференциальных уравнений при помощи разностных методов.
Изложение метода
При численном дифференцировании функцию аппроксимируют легко вычисляемой функцией и приближенно полагают . При этом можно использовать различные способы аппроксимации.
Интерполирование полиномами Ньютона
Рассмотрим случай аппроксимации интерполяционным многочленом Ньютона.
Пусть задана сетка - исследуемая функция. Введем обозначение и введем понятие разделенные разности , определяемые следующим образом:
- ,
- и т.д.
Можно доказать, что
- .
Запишем интерполяционный многочлен Ньютона и продифференцируем его почленно:
Общая формула примет следующий вид:
Обрывая ряд на некотором числе членов, получим приближенное выражение для соответствующей производной. Наиболее простые выражения получим, оставляя в формуле (1) только первый член:
- ,
- ,
Исследование точности полученных выражений при численных расчётах удобно делать при помощи апостериорной оценки, по скорости убывания членов ряда (1). Если шаг сетки достаточно мал, то погрешность близка к первому отброшенному члену. Пусть мы используем узлы . Тогда первый отброшенный член содержит разделенную разность , которая согласно (2) примерно равна . Перед ней стоит сумма произведений различных множителей ; каждое произведение содержит множителей, а вся сумма состоит из слагаемых. Отсюда следует оценка погрешности формулы (3) с узлами:
В частности, если сетка равномерная, то , откуда
- .
Стоит заметить, что строгое априорное исследование погрешности формулы (3), аналогичное выводу остаточного члена многочлена Ньютона в форме Коши, для произвольного расположения узлов приводит к той же оценке (3).
Таким образом, порядок точности формулы (1) по отношению к шагу сетки равен числу оставленных в ней членов, или, что то же самое, он равен числу узлов интерполяции минус порядок производной. Поэтому минимальное число узлов, необходимое для вычисления -й производной, равно ; оно приводит к формулам (2) и обеспечивает первый порядок точности. Эти выводы соответствуют общему принципу: при почленном дифференцировании ряда скорость его сходимости уменьшается.
Интерполирование полиномами Лагранжа
Рассмотрим неравномерную сетку и обозначим за , шаги этой сетки. В качества примера получим формулы численного дифференцирования, основанные на использовании многочлена Лагранжа , построенного для функции по трем точкам . Многочлен имеет вид
Отсюда получим
Это выражение можно принять за приближенное значение в любой точке ∈ . Его удобнее записать в виде , где , .
В частности, при получим , И если сетка равномерна, , то приходим к центральной разностной производной, . При использовании интерполяционного многочлена первой степени точно таким образом можно получить односторонние разностные производные и . Далее вычисляя вторую производную многочлена , получим приближенное выражение для при ∈:
≈
На равномерной сетке это выражение совпадает со второй разностной производной . Ясно, что для приближенного вычисления дальнейших производных уже недостаточно многочлена , надо привлекать многочлены более высокого порядка и тем самым увеличивать число узлов, участвующих в аппроксимации.
Порядок погрешности аппроксимации зависит как от порядка интерполяционного многочлена, так и от расположения узлов интерполирования. Получим выражение для погрешности аппроксимации, возникающей при замене выражением . Будем считать, что ∈ и что величины имеют один и тот же порядок малости при измельчении сетки. По формуле Тейлора в предположении ограниченности получим ,
где ,±
Отсюда приходим к следующим разложениям разностных отношений
Подставляя полученные формулы в выражение для разностной производной и приводя подобные слагаемые получим
, ∈ .
Отсюда видно,что разностное выражение аппроксимирует со вторым порядком.
Если подставить полученные ранее разностные отношения в выражение для второй производной многочлена , то имеем
Из этого выражения видно, что даже на равномерной сетке,т.е. когда , второй порядок аппроксимации имеет место лишь в точке , а относительно других точек (например,) выполняется аппроксимация только первого порядка. Таким образом, получим аппроксимацию лишь первого порядка.
Интеполирование кубическими сплайнами
Для того, чтобы избежать больших погрешностей в процессе приближения, весь отрезок [a,b] разбивают на частичные отрезки и на каждом из частичных отрезков приближенно заменяют функцию многочленом невысокой степени (так называемая кусочно-полиномиальная интерполяция). Одним из способов интерполирования на всем отрезке является интерполирование с помощью сплайн-функций.
Сплайн-функцией или сплайном называют кусочно-полиномиальную функцию, определенную на отрезке [a,b] и имеющую на этом отрезке некоторое число непрерывных производных. Преимущество сплайнов перед обычной интерполяцией является, во-первых их сходимость, во-вторых, устойчивость процесса вычислений.
Построение кубического сплайна.
Пусть на [a,b] задана непрерывная функция . Введем сетку и обозначим , . Сплайном соответствующим данной функции и данным узлам называется функция , удовлетворяющая следующим условиям:
а) на каждом сегменте , функция является многочленом третьей степени;
б) функция , а также её первая и вторая производная производные непрерывны на [a,b];
в) ,;
На каждом из отрезков , будем искать функцию в виде многочлена третьей степени
,
где , , где - коэффициенты, подлежащие определению. Доказано, что существует единственный кубический сплайн, определяемый условиями а)-в) и граничными условиями .
Для их нахождения используются следующие формулы
1) ,
2) Для определения коэффициентов получаем систему уравнений
,
(система решается методом прогонки) По найденным коэффициентам коэффициенты , определяются с помощью явных формул
3)
4)
Найдем производные введенного кубического сплайна, имеем
Рассмотрим оценку погрешности метода, которая зависит от выбора сеток и от гладкости . Для простоты изложения допустим, что сетка равномерная, т.е.
с шагом
От функции будем требовать существования непрерывной на [a,b] четвертой производной, ∈ . Кроме того, предположим, что выполнены граничные условия и такие же условия для сплайнов. Обозначим,
,
Пусть - кубический сплайн, построенный для функции на сетке . В следующей теореме приведены оценки погрешности интерполяции для функции и её производных ,
Теорема
Для ∈ справедливы оценки
Из этих оценок следует, что при (т.е. при ) последовательности , сходятся соответственно к функциям .
Обычно дифференцирование кубического сплайна позволят определить первую и вторую производную интерполяционного многочлена с хорошей точностью. Если надо вычислить более высокие производные, то целесообразно строить сплайны высоких порядков. Из-за большей трудоемкости этот способ редко используется. Способ дифференцирования с помощью сплайновой интерполяцией теоретически мало исследован.
Тригонометрическая интерполяция
Не всякую функцию целесообразно приближать алгебраическими многочленами. Рассмотрим тригонометрическую интерполяцию. Если - периодическая функция с периодом l, то естественно строить приближения с помощью функций
Таким образом, тригонометрическая интерполяция состоит в замене тригонометрическим многочленом
- ,
коэффициенты которого отыскиваются из системы уравнений
,
где .
Рекомендации программисту
При реализации достаточно сложно и трудоемко использовать методы сплайнов и метод тригонометрической интерполяции. Поэтому был рассмотрен случай интерполирования полиномом Лагранжа по трем точкам. При вычислении полиномов ошибка была ничтожна. При попытке в качестве функции взять синус,ошибка достигала очень большх цифр и был порядка 1.0, что достаточно велико для этой функции. Поэтому следует внимательно выбирать шаг, расстояние между точками. Программа требует ввода координат трех точек для полинома в порядке возрастания,а также точки промежуточной, в которой будет находиться значение производной. На экран выводится значение, полученное при помощи полинома Лагранжа,истинное значение, полученное аналитически и ошибка вычисления.
Список литературы
- А.А.Самарский, А.В.Гулин. Численные методы. Москва «Наука», 1989.
- Н.Н.Калиткин. Численные методы. Москва «Наука», 1978.