Интерполяционные многочлены в форме Лежандра и Ньютона

Интерполяционные многочлены в форме Лежандра и Ньютона Интерполяционный многочлен Лагранжа Интерполяционный многочлен Лагранжа имеет вид (39) или в компактной форме В случае равноотстоящих узлов xk=x0+k×h значение x можно представить в виде x = x0 +  t× h и многочлен Лагранжа записать как Представления (39)­(40) удобны при одиночных вычислениях, массовое же их использование достаточно трудоемко, кроме линейной (40) и квадратичной интерполяции (41) (42) Можно показать, что для функций, имеющих в данном диапазоне непрерывные производные  до (n+1)­го порядка, остаточный член интерполяционного многочлена (погрешность  аппроксимации) не превышает Rn(f)=  (43) Интерполяционный многочлен Лагранжа в случае монотонной функции может быть  использован и для решения обратной задачи интерполяции: поиска значения аргумента для заданного значения функции (44)В случае немонотонной функции, для которой обратная функция неоднозначна, приходится  строить интерполяционный многочлен и решать уравнение Ln(x)=f в выбранном диапазоне. 7.6.2.Конечные разности Пусть имеется таблица значений fk функции F(x) для равноотстоящих значений  аргумента xk = x0+k×h ( k=0, 1, 2,..,N ). Величины D fk= fk+1 ­ fk ( k=0, 1, 2,..,N­1 ) называются конечными разностями первого порядка, D2fk =D fk+1 ­D fk = fk+2 ­2 fk+1 +fk (k=0, 1, 2,..,N­2) ­ конечными разностями второго порядка и т.д. Обычно конечные разности записывают в виде таблиц; например, (45) (46) xk fk Df k D2fk D3fk D4fk 0 0 0 1 2 4 1 3 6 4 2 9 10 4 3 19 14 4 33       0         Можно показать, что отношение Dmfk / hk может быть принято за оценку m­й производной функции в точке xk . Возьмем для примера m=2. Разложим f(xk+2h) и f(xk+h) в ряд Тейлора в окрестности xk и, подставив в D2fk , получаем D2fk=fk+2 ­2 fk+1 +fk = h2 f ¢¢(xk) + h3 f¢¢¢(xk) + ... ,откуда xk fk Df k D2fk 0 1 2 4 1 3 6 4 2 9 10 4 3 19 14 4 4 33 18 4 5 51 22   6 73     f¢¢(xk) =  где O(h) ­ величина порядка h. Очевидно, что для функции ­ многочлена m­го порядка конечные разности m ­го порядка  должны быть одинаковыми, а m+1 ­ е и более высоких порядков ­ нулевыми. Эти  соображения могут быть положены в основу табличной экстраполяции (расширения). Так,  обнаружив в приведенной таблице постоянство вторых разностей обратным ходом  пополняем ее за пределами исходного диапазона. e ­5 e 15 e 0 e ­3  e 3 e ­4  e 6 e ­4  e 10 e ­10  e 5 e 0 0 0 0 0 0 0 0 e ­2  e 0 e e ­  e e ­ e e 0 0 0 0       ­20  e 15 e0 0 0 0 0   0                             При практическом использовании конечных разностей следует учитывать быстрый рост  погрешности. Пусть в одном из значений функции допущена погрешность порядка e.  Приведенная здесь таблица показывает, что конечная разность 6 ­го порядка содержит  ошибку, в 20 раз превосходящую исходную. Поэтому пользоваться разностями высших  порядков нужно весьма осторожно. 7.6.3.Интерполяционные формулы (равномерная сетка) Интерполяционные формулы (вывод их можно обнаружить практически во всех  руководствах по численному анализу) позволяют отыскивать значения табличной функции в точках, отличных от узлов таблицы, без построения интерполяционного многочлена. Самыми популярными из них являются интерполяционные формулы Ньютона. Формула Ньютона интерполирования вперед Формула Ньютона интерполирования назад Пример. Пусть задана функция и ее конечные разности [3] x k f k Df k D2f k D3f k D4f k 0.1 0.09983 0.09884 ­0.00199 ­0.00096 0.00002 0.2 0.19867 0.09685 ­0.00295 ­0.00094   0.3 0.29552 0.09390 ­0.00389   0.4 0.38942 0.09001   0.5 0.47943При поиске f(0,14) разумнее выбрать начальный узел x=0.1 (h=0.4) и воспользоваться  интерполированием вперед При поиске f(0,45) разумнее выбрать начальный узел x=0.5 (h=0.5) и воспользоваться  интерполированием назад При интерполировании в середине таблиц можно пользоваться и другими  интерполяционными формулами, которые строятся на основе конечных разностей,  последовательно выбираемых из выделенных клеток приведенной таблицы.   xk­3 f k­3 Df k­3 D2f k­3 D3f k­3 D4f k­3 D5f k­3   x k­2 f k­2 Dfk­2 D2f k­2 D3f k­2 D4f k­2 D5fk­2   xk­1 f k­1 Df k­1 D2f k­1 D3f k­1 D4f k­1     xk fk Dfk D2fk D3fk   xk+1 f k+1 Df k+1 D2f k+1   xk+2 f k+2 Df k+2   x k+3 fk+3                 Примерами таких формул могут служить следующие представления, где 0 < t < 1 : интерполяционные формулы Гаусса интерполяционная формула Стирлинга интерполяционная формула БесселяСуществуют и другие представления этих формул с использованием т.н. центральных  разностей, первая из которых определяется формулой d2f(x)= f(x+h/2) ­ f(x­h/2) . 7.6.4.Интерполирование на неравномерной сетке Выше мы утверждали, что погрешность интерполяционного многочлена Лагранжа  определяется величиной Rn(f)=  Если максимум n+1­й производной нам неподвластен, то произведение значений (x­xi) в ряде случаев (мы знаем функцию и желаем заменить ее более простой) существенно зависит от  выбора узлов интерполирования. П.Л.Чебышев доказал, что для tÎ[­1,1] минимум Q(t достигается при выборе узлов – 1£ x0m – число узлов в итоговом массиве Fu значений; диапазон вариации аргумента для обоих массивов  одинаков) Функция spline(X,Y,Z) выполняет интерполяцию Y=Y(X) кубическим сплайном и вывод  соответствующих значений в точках Z. Для получения большей информации используется  конструкция pp=spline(X,Y): здесь последующей командой V=ppval(pp,Z) можно найти  значения в точках Z, a командой [Xs, Coef, m,L]=unmkpp(pp) получить данные о векторе  разбиений аргумента Xs и его размерности m+1, о коэффициентах Coef, количество  которых равно mL (L=4). Функции interp1(X,Y,Z), interp1(X,Y,Z,’method’) обеспечивают одномерную табличную  интерполяцию для Y=Y(X) на массиве значений Z (если Y содержит несколько столбцов,  интерполяция ведется по каждому столбцу; значения Z должны входить в диапазон значений Х). Можно указать метод интерполяции – кусочно­линейной (linear, по умолчанию),  ступенчатой (nearest), кубической (cubic), кубическими сплайнами (spline).  Функцияinterp1q(X,Y,Z) реализует быструю линейную интерполяцию на неравномерной  сетке. Имеется возможность простой реализации двумерной табличной  интерполяции Y=Y(X1,X2) функциями interp2(X1,X2,Y,Z1,Z2),  interp2(X1,X2,Y,Z1,Z2,’method’) или griddata(X1,X2,Y,Z1,Z2), griddata(X1,X2,Y,Z1,Z2,  ’method’) для неравномерной сетки, трехмерной  интерполяцииY=Y(X1,X2,X3) – interp3(X1,X2,X3,Y,Z1,Z2, Z3), interp3(...,  ’method’), многомерной интерполяции Y=Y(X1,X2,...) ­ interpn (X1,X2,..,Y, Z1,Z2,..),  interp3(..., ’method’). Значения аргумента для исходных таблиц должны меняться  монотонно и задаваться в специальном формате функции meshgrid. Например, при желании построить график функции заданной при x1 и x2 в диапазонах от –1 до 1 на сравнительно крупной сетке, можно  выполнить интерполяцию на мелкой сетке и нарисовать график. Так, выполнив командыРис.7 » [X1,X2]=meshgrid(­1:0.1:1); % Исходная сетка » Y=exp(­X1.^2­X2.^2).*(1+X1+X2); % Значения функции » [Z1,Z2]=meshgrid(­1:0.05:1); % Сетка для интерполяции » Y2=interp2(X1,X2,Y,Z1,Z2); % Итог интерполяции » mesh(X1,X2,Y),hold on,mesh(Z1,Z2,Y2+2) % Графика получаем рис.7 (график функции на исходной сетке и смещенный на 2 вверх график итогов  интерполяции). Разумеется, для столь простой функции можно было не прибегать к  интерполяции.