Контрольная работа по "Математике"

• Запрос цены
• Пробная версия
• Запрос демонстрации

Интерполяция  функций интерполяционными  полиномами

В вычислительной математике существенную роль играет интерполяция функций, т.е. построение по заданной функции другой (как  правило, более простой), значения которой  совпадают со значениями заданной функции  в некотором числе точек. Причем интерполяция имеет как практическое, так и теоретическое значение. На практике часто возникает задача о восстановлении непрерывной функции  по ее табличным значениям, например полученным в ходе некоторого эксперимента. Для вычисления многих функций оказывается эффективно приблизить их полиномами или дробно-рациональными функциями (см., например [1] Ссылка в списке литературы). Теория интерполирования используется при построении и исследовании квадратурных формул для численного интегрирования, для получения методов решения дифференциальных и интегральных уравнений.

Все перечисленные  выше вопросы рассмотрены в классических учебниках по численным методам (см., например, [2-5] Ссылки в списке литературы). Цель этого раздела - демонстрация возможностей MATLAB для изучения вопросов, возникающих  при интерполяции функций, в основном при помощи интерполяционных полиномов. В данном разделе приводятся необходимые  сведения об интерполяции функций и  при помощи небольших программ, написанных на языке пакета MATLAB, изучаются проблемы, возникающие при интерполяции функций. Простота языка пакета MATLAB в сочетании  с широким набором его функций, в том числе и графических, позволяет вместо написания собственных  программ интерполяции и визуализации результатов сосредоточиться на исследовании большого числа примеров, что может быть использовано при  проведении лабораторных работ по численным  методам для студентов технических  факультетов вузов и институтов.

Задача  интерполяции функции, интерполяционные полиномы

Пусть на отрезке [a,b] задана функция ƒ(x). Задача интерполяции (или интерполирования) состоит в построении функции g(x), совпадающей с заданной ƒ(x) в некотором наборе точек {x₁,x₂,...,x_n+1} из отрезка [a,b] (эти точки называются узлами интерполяции), т.е. должны выполняться условия:

g(x_k)=y_k, k=1,2,...,n+1,

где y_k - известные значения функции ƒ(x) в точках x_k. Функция g(x) называется интерполянтом функции ƒ(x).

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

из которого видно, что узлы интерполяции не обязательно должны располагаться равномерно на отрезке [a,b].

Если ƒ(x) табличная функция, скажем полученная из эксперимента, т.е. известны только ее значения y_k в точках x_k, то, вообще говоря, о качестве полученного приближения судить трудно. Однако, если значения ƒ(x) могут быть вычислены в любой точке отрезка [a,b], то в этом случае можно исследовать качество получающегося приближения, например найдя максимальное уклонение функции g(x) от функции ƒ(x). На качество приближения сильное влияние оказывает количество и расположение узлов, а также гладкость функции ƒ(x). Эти вопросы и будут численно исследованы в следующих разделах.

Мы рассмотрим только линейную интерполяцию, т.е. такую, при которой функция g(x) разыскивается в виде линейной комбинации некоторых функций

где для k=1,2,...,n+1: φ_k(x) - заданные функции, а a_k - искомые коэффициенты.

Ясно, что из постановки задачи интерполяции (т.е. из совпадения значений интерполянта g(x) и интерполируемой функции ƒ(x) в точках x_k) следует, что коэффициенты a_k определяются из решения следующей системы линейных алгебраических уравнений:

или в развернутой  форме

Совершенно ясно, почему число коэффициентов a_k должно совпадать с числом узлов интерполяции x_k. Это нужно для того, чтобы матрица системы была квадратной (т.е. число неизвестных совпадало бы с числом условий, из которых находятся эти неизвестные). Кроме того, для однозначной разрешимости данной системы (при произвольной правой части) необходимо и достаточно, чтобы ее определитель был отличен от нуля, т.е.

Очень часто  в качестве системы функций  φ_k(x) выбирают полиномы, например степени x, именно:

Тогда соответствующий интерполянт называют интерполяционным полиномом. Существование и единственность интерполяционного полинома гарантируется, если все узлы интерполяции x_k различны. Действительно, определитель системы линейных алгебраических уравнений для нахождения коэффициентов a_k

является определителем  Вандермонда, который, как известно, равен

и, следовательно, отличен от нуля в случае, когда все узлы x_k различны. Поскольку матрица системы невырождена, то решение системы существует и единственно. Итак, интерполяционный полином существует и единственный.

Можно рассуждать и по-другому. Предположим, что есть два интерполяционных полинома g_k и h_k степени n такие, что для произвольного набора значений выполняются равенства g(x_k) = h(x_k) = y_k) для всех k=1,2,...,n+1, т.е. для n+1-ой точки. Тогда их разность h_k - g_k является полиномом степени не выше n, но обращается в ноль в n+1-ой точке. По известной теореме алгебры у полинома степени n не может быть больше чем n корней, следовательно h_k - g_k ≡ 0 и h_k ≡ g_k. Единственность установлена. А так как полином единственный, то у соответствующей системы линейных алгебраических уравнений есть только одно решение (для произвольной правой части). Из результатов линейной алгебры известно, что у системы может быть либо бесконечное число решений при некоторых правых частях, либо единственное, для произвольной правой части. Последнее как раз и выполняется.

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

Казалось бы, при практической реализации интерполяционного  процесса коэффициенты интерполяционного  полинома a_k можно найти, непосредственно решая систему линейных алгебраических уравнений

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

Установим зависимость  числа обусловленности этой матрицы  от числа узлов интерполяции, предполагая, что они распределены равномерно на отрезке [-1,1]. Для этого напишем  небольшую программу в MATLAB, в которой  в цикле по числу узлов генерируются матрицы, находится их число обусловленности  при помощи функции cond и строится график зависимости числа обусловленности матрицы от количества узлов

% задание границ  отрезка

a=-1; b=1;

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

N=1:30;

% инициализируем  массив (пока пустой) для записи  чисел обусловленности

C=[];

% проходим по  узлам в цикле

for n=N

    % формируем матрицу

    x=linspace(a,b,n)';

    M=[];

    for k=0:n

        M=[M x.^k];

    end

    % считаем  число обусловленности для текущего  числа узлов и записываем в  массив C

    C=[C cond(M)];

end

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

% в полулогарифмическом  масштабе

figure

semilogy(N,C)

grid on

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

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

Важно, что какие  бы подходы для построения интерполяционного  полинома не применялись, они всегда должны привести к одинаковому результату (если все вычисления проводятся точно, а не на компьютере), поскольку интерполяционный полином степени существует и  единственный при различных -ом узлах интерполяции. Другое дело, что разные способы построения интерполяционного полинома могут обладать разными вычислительными свойствами. Рассмотрим сначала интерполяционный полином в форме Лагранжа. Далее мы будем использовать обозначение для интерполяционного полинома в зависимости от способа его построения.

Интерполяционный  полином в форме  Лагранжа

Вывод формулы

Итак, мы ищем полином  L_n(x) степени не выше n, значения которого совпадают со значениями y_k заданной функции ƒ(x) в узлах x_k, где k=1,2,...,n+1 и все узлы различны.

Одним из способов записи интерполяционного полинома является форма Лагранжа. Предположим, что для k=1,2,...,n+1 функции Ф_n(x) являются полиномами степени n, которые обладают следующим свойством

Тогда полином

будет как раз  тем, который нам и нужен, поскольку  это полином степени не выше n и L_n(x_k) = y_k для всех k=1,2,...,n+1.

Функции Ф_n(x) строятся легко. Действительно, функция

является полиномом  степени n, который обращается в ноль для всех x_j не равных x_k. В точке x_k она принимает значение

Тогда

или, что то же самое:

Для примера  построим графики функций Ф_n(x) для пяти узлов, равномерно распределенных на отрезке [0,1], для чего можно воспользоваться следующей небольшой программой:

% задаем границы  отрезка и значение n

a=0; b=1; n=4

% генерируем набор  из n+1 точек, равномерно распределенных  на [0,1]

x=linspace(a,b,n+1);

% генерируем вспомогательный  набор из 1000 точек, равномерно  распределенных на [0,1]

% для вычисления  в них функций Psi_k для построения их графиков

xx=linspace(a,b,1000);

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

figure

axes('NextPlot','add')

grid on

% в цикле вычисляем  значения каждой функции Psi_k в точках, записанных в xx

% и строим графики  функции Psi_k

for k=1:n+1

    Phi=ones(size(xx));

    J=[1:k-1 k+1:n+1];

    for j=J

        Phi=Phi.*(xx-x(j))/(x(k)-x(j));

    end

    plot(xx,Phi,'Color',[rand rand rand],'LineWidth',2)

end

% отображаем узлы  интерполяции красными маркерами

plot(x,zeros(size(x)),'LineStyle','none','Marker','.','Color','r','MarkerSize',20)

На следующем  рисунке приведены полученные графики  Ф_n(x) для k=1,2,3,4,5, узлы отмечены красными маркерами. Видно, что каждая из функций Ф_n(x) равна единице только в одном из узлов интерполяции, а в остальных она равна нулю.

Итак, интерполяционный полином в форме Лагранжа имеет  вид

или в развернутой  форме:

Вычисление  интерполяционного  полинома в форме  Лагранжа

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

то для вычисления каждого из n+1-го произведений

требуется n делений, n-1 умножение и 2n сложений (вычитаний). Далее, когда произведения найдены, для вычисления суммы требуется n+1 умножений и n сложений. Т.е. получается, что для каждого значения x вычисление L_n(x) требует порядка n² операций.

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

за n-1 умножений, n сложений и 1 деление (эта операция делается для данного интерполяционного полинома только один раз перед вычислением его значений в точках), а затем находить значения интерполяционного полинома по формуле

L_n(x) = w(x)s(x),

где

При таком способе  вычисления значения интерполяционного  полинома в точке требуется предварительная  работа порядка n² операций, зато для каждого значения x время вычисления интерполяционного полинома требует порядка n операций. В знаменателе в выражении для функции s(x) находится разность x - x_k, однако при x = x_k значение интерполяционного полинома вычислять не требуется, оно известно и равно y_k.

Для вычисления интерполяционного полинома в MATLAB можно  применить оба вышеописанных  способа. Ниже приведены тексты функций, в которых реализованы эти  способы:

• в функции lagrange - первый способ;
• в функции lagrangef - второй способ;

Входные массивы x и y должны содержать значения x_k и y_k, соответственно, для всех k=1,2,...,n+1. Входной аргумент xx функций lagrange и lagrangef может быть массивом значений аргумента, для которых требуется вычислить значение интерполяционного полинома. Тогда в выходном аргументе yy вернется массив соответствующих значений полинома. При программировании функций lagrange и lagrangef не потребовалось делать цикла по элементам массива xx благодаря поддержке поэлементных операции при работе с массивами в MATLAB.

Примечание (про  поэлементные операции). 
Для массивов в MATLAB допустимы поэлементные умножение, деление и возведение в степень. Для этих операций применяются, соответственно, следующие сочетания (с точкой): .* , ./ , .^. Приведем несколько примеров поэлементных операций с векторами:

>> a=[10 20 30];

>> b=[2 4 6];

>> c=a./b

c =

     5     5     5 

>> d=a.*b

d =

    20    80   180

При программировании интерполяционного полинома по второму  способу в функции lagrangef не делалась проверка на равенство x какому-либо узлу, поскольку в MATLAB операция деления на ноль допустима (при делении на ноль числа не равного нулю получается Inf, а при делении нуля на ноль получается NaN, т.е Not a Number, не число).

function yy=lagrange(x,y,xx)

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

% x - массив координат  узлов

% y - массив значений  интерполируемой функции

% xx - массив значений аргумента, для которых надо вычислить значения полинома

% yy - массив значений полинома в точках xx  

% узнаем число  узлов интерполяции (N=n+1)

N=length(x);

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

yy=zeros(size(xx));

% в цикле считаем  сумму по узлам

for k=1:N

    % вычисляем  произведения, т.е. функции Psi_k

    t=ones(size(xx));

    for j=[1:k-1, k+1:N]

        t=t.*(xx-x(j))/(x(k)-x(j));

    end

    % накапливаем  сумму

    yy = yy + y(k)*t;

end 

function yy=lagrangef(x,y,xx)

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

% x - массив координат  узлов

% y - массив значений  интерполируемой функции

% xx - массив значений аргумента, для которых надо вычислить значения полинома

% yy - массив значений полинома в точках xx  

% узнаем число  узлов интерполяции (N=n+1)

N=length(x);

% предварительное  вычисление значений z_k

z=zeros(size(x));

for k=1:N

    t=1;

    for j=[1:k-1 k+1:N]

        t=t*(x(k)-x(j));

    end

    z(k)=y(k)/t;

end

% вычисление w(x)

w=ones(size(xx));

for k=1:N

    w=w.*(xx-x(k));

end

s=zeros(size(xx));

% вычисление s(x)

for k=1:N

    s=s+z(k)./(xx-x(k));

end

% вычисление значений  интерполяционного полинома

yy=w.*s;

В качестве тестового  примера проинтерполируем функцию sin x на отрезке [1,9] с шагом 2 и построим графики sin x и полученного интерполяционного полинома. В данном случае будет 5 узлов и, следовательно, интерполяционный полином получится 4-ой степени. Для построения графика интерполяционного полинома вычислим его значения при помощи функции lagrange в 1000 равномерно отстоящих друг от друга точках на отрезке [1,9].

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

x=1:2:9;

y=sin(x);

% задание точек,  в которых требуется найти  значения интерполяционного полинома

xx=linspace(1,9,1000);

% нахождение значений  интерполяционного полинома

yy=lagrange(x,y,xx);

% построение графиков

figure('Color','w')

% вывод графика  sin(x)

fplot(@sin,[1 9])

hold on

% вывод графика  полинома

plot(xx,yy,'r')

% вывод узлов  интерполяции

plot(x,y,'bo')

% размещение легенды

legend('sin\itx','{\itL_n}({\itx})','nodes',-1)

Получающийся  результат приведен на рисунке ниже

Хорошим тестовым примером для функций lagrange и lagrangef является приближение полиномиальной функции. Выберем в качестве интерполируемой функции полином пятой степени

p₅(x) = 4x⁵ - 3x⁴ + 14x³ - 22x² - x + 5

и проинтерполируем его полиномом пятой степени L₅(x) на отрезке [-1, 1.5]. Для интерполяции понадобится шесть узлов, мы выберем их равноотстоящими на этом отрезке и проведем интерполяцию при помощи интерполяционного полинома в форме Лагранжа, вычисляемого приведенной выше функцией lagrange. Мы выведем не только график исходного полинома p₅(x) и график интерполяционного полинома L₅(x), но и график ошибки, т.е. разности p₅(x) - L₅(x).

% выбор 6 равноотстоящих  узлов на отрезке [-1, 1.5]

x=linspace(-1,1.5,6);

% задание inline-функции для вычисления полинома

p5=inline('4*x.^5 - 3*x.^4 + 14*x.^3 - 22*x.^2 - x + 5');

% вычисление полинома p5 в узлах интерполяции

y=p5(x);

% задание 10000 равноотстоящих  точек на отрезке [-1, 1.5] для  вычисления в них значений

% интерполяционного  полинома

xx=linspace(-1,1.5,10000);

% вычисление значений  интерполяционного полинома

yy=lagrange(x,y,xx);

% графический вывод  результатов 

figure('Color','w')

subplot(2,1,1)

% вывод графика  исходного полинома p5

fplot(p5,[-1 1.5])

hold on

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

plot(xx,yy,'r')

% вывод легенды

legend('{\itp}_5({\itx})','{\itL}_5({\itx})',-1)

subplot(2,1,2)

% вывод графика  распределения ошибки на отрезке  [-1, 1.5]

plot(xx,p(xx)-yy)

% вывод заголовка  к графику для ошибки

title('error = {\itp}_5({\itx}) - {\itL}_5({\itx})')

В результате получаем, что график интерполяционного полинома пятой степени практически совпадает  с графиком исходного полинома, однако ошибка вовсе не равна нулю. Она  имеет порядок 10^-14, что объясняется ошибками округлений при вычислениях значений интерполяционного полинома. Более того, график распределения ошибки имеет характерный вид - в середине интервала ошибка меньше, чем по краям.

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

Проведем еще  ряд экспериментов. Приблизим полином  пятой степени интерполяционным полиномом с числом узлов меньше шести и с числом узлов, значительно большим, чем 6.

Сначала проинтерполируем полином

p₅(x) = 4x⁵ - 3x⁴ + 14x³ - 22x² - x + 5

на отрезке [-1, 1.5] с тремя равноотстоящими узлами и получим интерполяционный полином  L₂(x) второй степени, затем с четырьмя и пятью равноотстоящими узлами и получим, соответственно, полиномы L₃(x) и L₄(x) третьей и четвертой степени.

% задание inline-функции для вычисления полинома

p5=inline('4*x.^5 - 3*x.^4 + 14*x.^3 - 22*x.^2 - x + 5');

% задание трех  равноотстоящих узлов на [-1, 1.5]

x3=linspace(-1,1.5,3);

% вычисление в  них исходного полинома p5(x)

y3=p5(x3);

% задание четырех  равноотстоящих узлов на [-1, 1.5]

x4=linspace(-1,1.5,4);

% вычисление в  них исходного полинома p5(x)

y4=p5(x4);

% задание пяти  равноотстоящих узлов на [-1, 1.5]

x5=linspace(-1,1.5,5);

% вычисление в  них исходного полинома p5(x)

y5=p5(x5);

% задание 10000 узлов  на [-1, 1.5] для вычисления в них  значений интерполяционных полиномов

xx=linspace(-1,1.5,10000);

% вычисление в  них значений интерполяционных  полиномов 2-ой, 3-ей и 4-ой степени

yy3=lagrange(x3,y3,xx);

yy4=lagrange(x4,y4,xx);

yy5=lagrange(x5,y5,xx);

% графический вывод  результатов

figure('Color','w')

% построение графика  исходного полинома p5(x)

fplot(p5,[-1 1.5],'r')

hold on

% построение графиков  интерполяционных полиномов 2-ой, 3-ей и 4-ой степени

plot(xx,yy3,'k')

plot(xx,yy4,'b')

plot(xx,yy5,'g')

% вывод легенды

legend('{\itp}_5({\itx})','{\itL}_2({\itx})',...

    '{\itL}_3({\itx})','{\itL}_4({\itx})',-1)

В результате получаем вполне предсказуемый результат, чем  ближе степень интерполяционного  полинома к 5-ой, тем лучше приближение.

Теперь возьмем  интерполяционные полиномы достаточно высокой степени и проследим  за ошибкой интерполирования. Для  примера проинтерполируем полином

p₅(x) = 4x⁵ - 3x⁴ + 14x³ - 22x² - x + 5

на отрезке [-1, 1.5] с десятью, двадцатью и тридцатью  равноотстоящими узлами и получим, интерполяционные полиномы L₉(x), L₁₉(x) и L₂₉(x), соответственно, девятой, девятнадцатой и двадцать девятой степени.

% задание inline-функции для вычисления полинома

p5=inline('4*x.^5 - 3*x.^4 + 14*x.^3 - 22*x.^2 - x + 5');

% задание десяти  равноотстоящих узлов на [-1, 1.5]

x10=linspace(-1,1.5,10);

% вычисление в  них исходного полинома p5(x)

y10=p5(x10);

% задание двадцати  равноотстоящих узлов на [-1, 1.5]

x20=linspace(-1,1.5,20);

% вычисление в  них исходного полинома p5(x)

y20=p5(x20);

% задание тридцати  равноотстоящих узлов на [-1, 1.5]

x30=linspace(-1,1.5,30);

% вычисление в  них исходного полинома p5(x)

y30=p5(x30);

% задание 10000 узлов  на [-1, 1.5] для вычисления в них  значений интерполяционных полиномов

xx=linspace(-1,1.5,10000);

% вычисление в  них значений интерполяционных  полиномов 9-ой, 19-ой и 29-ой степени

yy10=lagrange(x10,y10,xx);

yy20=lagrange(x20,y20,xx);

yy30=lagrange(x30,y30,xx);

% графический вывод  результатов

figure('Color','w')

% вывод ошибки  интерполирования для интерполяционного  полинома 9-ой степени

subplot(3,1,1)

plot(xx,yy10-p5(xx),'b')

title('error = {\itp}_5({\itx}) - {\itL}_{9}({\itx})')

% вывод ошибки  интерполирования для интерполяционного  полинома 19-ой степени

subplot(3,1,2)

plot(xx,yy20-p5(xx),'g')

title('error = {\itp}_5({\itx}) - {\itL}_{19}({\itx})')

% вывод ошибки  интерполирования для интерполяционного  полинома 29-ой степени

subplot(3,1,3)

plot(xx,yy30-p5(xx),'r')

title('error = {\itp}_5({\itx}) - {\itL}_{29}({\itx})')

В результате мы видим, что ошибки округлений при  вычислении значений интерполяционного  полинома приводят к тому, что ошибка интерполирования растет (порядка 10^-14 для интерполяционного полинома 9-ой степени, порядка 10^-11 для интерполяционного полинома 19-ой степени и порядка 10^-9 для интерполяционного полинома 29-ой степени), хотя она должна была бы равняться нулю (при вычислениях в точной арифметике). Заметим, что характер распределения ошибки сохраняется - чем ближе к границам отрезка интерполирования, тем ошибка больше.

Вычислительные  свойства интерполяционного полинома делают его неприменимым при достаточно большом числе узлов интерполирования. Изучим зависимость максимального  уклонения интерполяционного полинома от интерполируемой функции снова  на примере полинома

p₅(x) = 4x⁵ - 3x⁴ + 14x³ - 22x² - x + 5

который мы будем интерполировать полиномами на отрезке [-1, 1.5] с числом узлов интерполирования, изменяющимся от двух до ста. Т.е. наша цель - получить график зависимости от степени интерполяционного полинома n. Для этого в цикле по числу узлов интерполяции будем строить интерполяционный полином и вычислять максимальное уклонение его от p₅(x), определяя это уклонение по 1000 равномерно отстоящим точкам на отрезке [-1, 1.5].

% задание inline-функции для вычисления полинома

p5=inline('4*x.^5 - 3*x.^4 + 14*x.^3 - 22*x.^2 - x + 5');

% инициализация  массива для записи ошибок  интерполяции 

err=[];

% задание массива  с количеством узлов интерполяции

K=2:100;

% построение в  цикле интерполяционных полиномов  и вычисление ошибки интерполяции

for k=K

    % задание  k равноотстоящих на отрезке [-1, 1.5] узлов интерполяции

    x=linspace(-1,1.5,k);

    % вычисление  в них значений интерполяционного  полинома

    y=p5(x);

    % задание  1000 равноотстоящих на отрезке  [-1, 1.5] точек

    xx=linspace(-1,1.5,1000);

    % нахождение  в них значений интерполяционного  полинома

    yy=lagrange(x,y,xx);

    % вычисление ошибки интерполяции

    e=max(abs(p5(xx)-yy));

    % добавление  ее в массив

    err=[err e];

end

% графический вывод  результатов

figure('Color','w')

% построение зависимости  ошибки от степени интерполяционного  полинома

semilogy(K-1,err,'LineWidth',2)

% нанесение сетки

grid on

% добавление заголовка

title('Dependence of the max|{\itL_n}({\itx})-{\itp}_5({\itx}))| on {\itn}')

В результате работы этой программы мы видим (см. рисунок  ниже), что сначала ошибка уменьшается  и достигает своего минимума (когда  степень интерполяционного полинома равна 5). Как уже обсуждалось выше, она не равна в точности нулю из-за ошибок округления при вычислении интерполяционного  полинома. Далее, при повышении степени  интерполяционного полинома, ошибка начинает быстро возрастать, что свидетельствует  о практической неприменимости интерполяционных полиномов высоких степеней.

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

Ошибка  интерполяции, чебышевские узлы

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

Погрешность интерполирования функции ƒ(x) интерполяционным полиномом L_n(x) n-ой степени в точке x выражается следующим образом [3]

где - некоторая точка из промежутка [a,b] и ω_n+1(x) - полином степени n + 1, определяемый по формуле

ω_n+1(x) = (x-x₁)(x-x₂)...(x-x_n+1)

где x₁,x₂,...,x_n+1 - узлы интерполяционного полинома L_n(x).

Значение ω вообще говоря неизвестно, однако, если известно максимальное значение на отрезке [a,b]:

(или его оценка  сверху), то можно пользоваться  следующей оценкой сверху для  ошибки интерполяции в точке  x:

Продемонстрируем  использование этой оценки для проверки ошибки интерполяции на примере интерполяции функции ƒ(x) = sinx на отрезке [-π,π] интерполяционным полиномом 4-ой степени с равномерно отстоящими пятью узлами интерполяции. Для вычисления интерполяционного полинома L₄(x) воспользуемся функцией lagrange, текст которой приведен в разделе Вычисление интерполяционного полинома в форме Лагранжа. Сначала зададим узлы интерполяции и вычислим в них значения функции sinx:

>> x=linspace(-pi,pi,5)

x =

   -3.1416   -1.5708         0    1.5708    3.1416

>> y=sin(x)

y =

   -0.0000   -1.0000         0    1.0000    0.0000

Теперь найдем значение yy интерполяционного полинома L₄(x) и ошибку интерполяции err в точке :

>> xx=pi/6

xx =

    0.5236

>> yy=lagrange(x,y,xx)

yy =

    0.4321

>> err=abs(sin(xx)-yy)

err =

    0.0679

Вычислим правую часть оценки est с учетом того, что производные функции sinx не превосходят единицу, т.е. заменим M_n+1 на 1:

>> est=abs(prod(x-xx))/factorial(5)

est =

    0.0918

Видим, что est (оценка) больше, чем err (ошибка интерполяции) при x = .

Мы проделали  проверку оценки ошибки в точке, а  следующие команды приводят к  построению графиков ошибки интерполяции функции ƒ(x) = sinx и ее оценки сверху на отрезке [-π,π] интерполяционным полиномом 4-ой степени с равномерно отстоящими пятью узлами интерполяции.

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

x=linspace(-pi,pi,5);

% вычисление в  них значений функции

y=sin(x);

% задание промежуточных  точек

xx=linspace(-pi,pi,300);

% вычисление в  них ошибки интерполяции

yy=lagrange(x,y,xx);

err=abs(yy-sin(xx));

% вычисление оценки  ошибки в промежуточных точках

est=zeros(size(xx));

for i=1:length(xx)

    est(i)=abs(prod(x-xx(i)))/factorial(5);

end

% построение графиков  ошибки и ее оценки сверху

figure('Color','w')

plot(xx,err,'g',xx,est,'r')

legend('error','estimate')

То, что в узлах  интерполяции ошибка интерполяции и ее оценка сверху обращаются в ноль, очевидно, т.к. значения интерполянта в узлах совпадают со значениями интерполируемой функции и правая часть оценки

обращается в  ноль, т.к.

ω_n(x) = (x-x₁)(x-x₂)...(x-x_n+1) ≡ 0

в узлах интерполяции x₁,x₂,...,x_n+1.

Если нас интересует максимальное уклонение интерполяционного  полинома от интерполируемой функции, то в последней оценке можно взять  максимум по x:

Из этой оценки видно, что мы можем управлять  максимальной ошибкой за счет подходящего  выбора узлов интерполяции. Цель выбора узлов интерполяции {x₁,x₂,...,x_n+1} - сделать максимальное значение модуля полинома

ω_n+1(x) = (x-x₁)(x-x₂)...(x-x_n+1) ≡ 0

на отрезке [a,b] как можно меньше. Разумеется, эта задача имеет смысл только если мы можем вычислять интерполируемую функцию ƒ(x) в любых точках на отрезке [a,b]. Если функция задана таблично, то такой возможности нет.

Ограничимся сначала  случаем a = -1, b = 1. Известно, что если узлы интерполяции {x₁,x₂,...,x_n+1} являются корнями полинома Чебышева степени n + 1, то величина принимает наименьшее возможное значение по сравнению с любым другим выбором набора узлов интерполяции.

Полиномы Чебышева, предложенные и исследованные П.Л. Чебышевым в середине 19-го века, определяются следующим образом:

T_k(x) = cos(k arccos x), |x|≤1.

Очевидно, что  для k = 1 функция T₁(x) действительно является полиномом первой степени, поскольку

T₁(x) = cos(arccos x) = x

В случае k = 2 тоже очевидно, что T₂(x) есть полином второй степени, если воспользоваться известным тригонометрическим тождеством

cos2θ = 2cos²θ - 1,

положив θ = arccos x. Тогда получаем

T₂(x) = 2x² - 1,

Тригонометрическое  тождество

cos(k + 1)θ = 2cosθcoskθ - cos(k - 1)

позволяет легко  установить, что для полиномов  Чебышева справедливо следующее  рекуррентное соотношение

T_k+1(x) = 2xT_k(x) - T_k-1(x)

При помощи этого  рекуррентного соотношения можно  последовательно получить формулы  для полиномов Чебышева любой  степени. Из этого рекуррентного  соотношения видно, что коэффициент  при старшей степени x^k равен 2^k-1.

В MATLAB полиномы Чебышева можно вывести при помощи Symbolic Math Toolbox, обратившись к функции Maple, которая называется ChebyshevT и функции simplify (для вызова функций Mаple из MATLAB предназначена функция maple):

>> p2=maple('simplify(ChebyshevT(2,x),''ChebyshevT'')')

p2 =

-1+2*x^2 

>> p3=maple('simplify(ChebyshevT(3,x),''ChebyshevT'')')

p3 =

4*x^3-3*x 

>> p4=maple('simplify(ChebyshevT(4,x),''ChebyshevT'')')

p4 =

1+8*x^4-8*x^2 

>> p5=maple('simplify(ChebyshevT(5,x),''ChebyshevT'')')

p5 =

16*x^5-20*x^3+5*x 

Построим  графики  полиномов Чебышева для нескольких различных значений  , для чего выполним следующие команды: 

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

figure('Color','w')

axes

grid on

hold on

% запись в строковую  переменную выражения для полинома  Чебышева 2-ой степени

p2=maple('simplify(ChebyshevT(2,x),''ChebyshevT'')');

% создание inline-функции

t2=inline(p2);

% построение ее  графика на отрезке [-1,1]

fplot(t2,[-1,1],'r')

% аналогично для  полиномов 4-ой, 7-ой и 10-степеней

p4=maple('simplify(ChebyshevT(4,x),''ChebyshevT'')');

t4=inline(p4);

fplot(t4,[-1,1],'g')

p7=maple('simplify(ChebyshevT(7,x),''ChebyshevT'')');

t7=inline(p7);

fplot(t7,[-1,1],'b')

p10=maple('simplify(ChebyshevT(10,x),''ChebyshevT'')');

t10=inline(p10);

fplot(t10,[-1,1],'m')

% вывод легенды

legend('{\itk}=2','{\itk}=4','{\itk}=7','{\itk}=10')

Корни полинома Чебышева легко найти, решив уравнение

из которого видно, что полином T_k(x) имеет ровно k различных корней, расположенных на отрезке [-1,1]:

которые и следует выбирать в качестве узлов интерполирования. Корни полиномов Чебышева расположены симметрично относительно нуля на отрезке [-1,1] и неравномерно - чем ближе к краям отрезка, тем корни расположены плотнее.

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

Если мы положим

для того, чтобы коэффициент при старшей степени полинома ω_k(x) был равен 1 (по определению полинома ω_k(x)), то получим, что

Несложно показать (ссылка в списке литературы), что  для любого полинома p_k(x) степени k с коэффициентом при старшей производной равным единице будет верно следующее неравенство

т.е. полиномы Чебышева являются полиномами наименее уклоняющимися от нуля.

Таким образом, выбор в качестве узлов интерполирования корней полинома Чебышева является наилучшим  в смысле минимизации правой части  оценки

которая теперь приобретает вид:

где

Продемонстрируем  преимущество выбора корней полиномов  Чебышева в качестве узлов интерполяции. Для примера будем интерполировать  функцию ƒ(x) = sinx на отрезке [-1,1] полиномами различных степеней для равноотстоящих и чебышевских узлов и сравнивать получающуюся ошибку интерполяции, строя графики ошибки L_n(x) - sin(x). Для этого можно воспользоваться последовательностью команд, приведенной ниже, в которой следует изменять значение k (т.е. число узлов интерполяции).

% задание числа узлов интерполяции

k=7;

% генерация равномерно  отстоящих узлов на отрезке  [-1,1]

x=linspace(-1,1,k);

% вычисление в  них значений функции синус

y=sin(x);

% генерация 500 равномерно  отстоящих точек на отрезке  [-1,1]

xx=linspace(-1,1,500);

% вычисление в  них значений интерполяционного  полинома

yy=lagrange(x,y,xx);

% нахождение ошибки  интерполяции в этих точках

err=yy-sin(xx); 

% вычисление корней  полинома Чебышева в качестве  узлов на отрезке [-1,1]

m=0:k-1;

xcheb=cos((2*m+1)/k*0.5*pi);

% вычисление в  них значений функции синус

ycheb=sin(xcheb);

% вычисление значений  интерполяционного полинома в  500 равномерно отстоящих 

% точек на отрезке  [-1,1]

yy=lagrange(xcheb,ycheb,xx);

% нахождение ошибки  интерполяции в этих точках

errcheb=yy-sin(xx);

% построение графиков  ошибок для равномерных и чебышевских узлов

figure('Color','w')

plot(xx,err,xx,errcheb)

% вывод легенды и заголовка

legend('equidistance','chebyshev')

title('Interpolation error')

Ниже приведены  графики ошибки интерполяции для  равноотстоящих и чебышевских узлов для нескольких чисел узлов интерполяции. На них отчетливо прослеживается типичная ситуация, когда при интерполяции с равноотстоящими узлами ошибка интерполяции больше на краях промежутка интерполяции, а при интерполяции с чебышевскими узлами ошибка более равномерно распределена на интервале интерполирования.

Мы рассмотрели  выбор чебышевских узлов, когда отрезок, на котором производится интерполирование функции, есть [-1,1]. Переход к общему случаю интерполирования функции на произвольном отрезке [a,b] не представляет никаких трудностей. Достаточно сделать линейную замену переменных, переводящую отрезок [-1,1] в отрезок [a,b]. Такая замена выражается следующей простой формулой:

Тогда для произвольного  отрезка [a,b] чебышевскими узлами интерполяции будут узлы, вычисляемые по формуле:

где m = 0,1,...,k-1. Здесь k - число узлов интерполяции, т.е. k = n + 1, где n - степень интерполяционного полинома L_n(x).

Оценка ошибки интерполяции для отрезка [a,b] приобретает вид:

где

Приведем пример интерполирования с чебышевскими и равноотстоящими узлами функции

на отрезке [-2,3]. Будем выводить график ошибки интерполирования для различного числа равноотстоящих и чебышевских узлов, для чего в приведенной ниже программе следует менять значение переменной k.

% задание границ  отрезка интерполирования

a=-2;

b=3;

% задание inline-функции

fun=inline('sin(x).^2+sin(x.^2)');

% задание числа  узлов интерполяции

k=7;

% генерация равномерно  отстоящих узлов на отрезке  [a,b]

x=linspace(a,b,k);

% вычисление в  них значений функции sin(x)^2+sin(x^2)

y=fun(x);

% генерация 500 равномерно  отстоящих точек на отрезке  [a,b]

xx=linspace(a,b,500);

% вычисление в  них значений интерполяционного  полинома

yy=lagrange(x,y,xx);

% нахождение ошибки  интерполяции в этих точках

err=yy-fun(xx); 

% вычисление чебышевских узлов на отрезке [a,b]

m=0:k-1;

xcheb=0.5*((b-a)*cos((2*m+1)/k*0.5*pi)+a+b);

% вычисление в  них значений функции sin(x)^2+sin(x^2)

ycheb=fun(xcheb);

% вычисление значений  интерполяционного полинома в  500 равномерно отстоящих 

% точек на отрезке  [a,b]

yy=lagrange(xcheb,ycheb,xx);

% нахождение ошибки  интерполяции в этих точках

errcheb=yy-fun(xx);

% построение графиков  ошибок для равномерных и чебышевских узлов

figure('Color','w')

plot(xx,err,xx,errcheb)

% вывод легенды и заголовка

legend('equidistance','chebyshev')

title('Interpolation error')

Ниже приведены  графики ошибки интерполяции для  равноотстоящих и чебышевских узлов для нескольких чисел узлов интерполяции. Ситуация аналогична случаю отрезка , рассмотренному выше. При интерполяции с равноотстоящими узлами ошибка интерполяции больше на краях промежутка интерполяции, а при интерполяции с чебышевскими узлами ошибка более равномерно распределена на интервале интерполирования и меньше по значению, чем при интерполяции с равноотстоящими узлами.

Сходимость  интерполяционного  полинома

Важно определить, что мы понимаем под сходимостью  интерполяционного полинома L_n(x) к интерполируемой функции ƒ(x).

Говорят, что  интерполяционный полином L_n(x) сходится к ƒ(x) в точке x, если

и интерполяционный полином L_n(x) равномерно сходится к ƒ(x), если

Отметим, что  из поточечной сходимости не следует  равномерная и, если нет поточечной сходимости, то нет и равномерной.

Как гласит теорема  Вейерштрасса, всякая непрерывная функция  на отрезке может быть равномерно приближена полиномом с любой  точностью, однако из этой теоремы не следует существование именно интерполяционного  полинома, который приближал бы с  любой точностью заданную непрерывную  функцию.

Мы начнем изучение сходимости интерполяционного полинома с классического примера Рунге (1901г.). Рунге исследовал интерполяцию полиномами функции

с равноотстоящими  узлами на отрезке [-1,1] и доказал, что  вблизи границ этого отрезка при  x>0.72 погрешность интерполяции ƒ(x)-L_n(x) неограниченно увеличивается с увеличением степени интерполяционного полинома n.

Проинтерполируем функцию Рунге r(x) на отрезке [-1,1], увеличивая число равноотстоящих узлов, и отобразим графически функцию r(x) и интерполяционные полиномы на одних осях, а погрешности интерполирования на других осях в одном графическом окне.

% задание функции  Рунге 

r=inline('1./(1+25*x.^2)');

% создание графического  окна и двух пар осей

figure('Color','w')

subplot(2,1,1)

fplot(r,[-1,1],'r')

hold on

subplot(2,1,2)

hold on

% генерация 500 равномерно  отстоящих точек на отрезке  [-1,1]

xx=linspace(-1,1,500); 

% задание 6 равноотстоящих  узлов интерполирования

k=6;

x=linspace(-1,1,k);

% вычисление в  них значения функции Рунге

y=r(x);

% вычисление значений  интерполяционного полинома в  точках xx

% и построение  его графика

yy=lagrange(x,y,xx);

subplot(2,1,1)

plot(xx,yy,'b')

% нахождение ошибки  интерполяции в этих точках

% и построение  ее графика

err=yy-r(xx);

subplot(2,1,2)

plot(xx,err,'b') 

% задание 10 равноотстоящих  узлов интерполирования

k=10;

x=linspace(-1,1,k);

% вычисление в  них значения функции Рунге

y=r(x);

% вычисление значений  интерполяционного полинома в  точках xx

% и построение  его графика

yy=lagrange(x,y,xx);

subplot(2,1,1)

plot(xx,yy,'g')

% нахождение ошибки  интерполяции в этих точках

% и построение  ее графика

err=yy-r(xx);

subplot(2,1,2)

plot(xx,err,'g') 

% задание 12 равноотстоящих  узлов интерполирования

k=12;

x=linspace(-1,1,k);

% вычисление в  них значения функции Рунге

y=r(x);

% вычисление значений  интерполяционного полинома в  точках xx

% и построение  его графика

yy=lagrange(x,y,xx);

subplot(2,1,1)

plot(xx,yy,'k')

% нахождение ошибки  интерполяции в этих точках

% и построение  ее графика

err=yy-r(xx);

subplot(2,1,2)

plot(xx,err,'k') 

% вывод легенд  на оси

subplot(2,1,1)

legend('{\itr}({\itx})', '{\itL}_6({\itx})', '{\itL}_{10}({\itx})', '{\itL}_{12}({\itx})',-1)

subplot(2,1,2)

legend('{\itL}_6({\itx})-{\itr}({\itx})', '{\itL}_{10}({\itx})-{\itr}({\itx})',...

    '{\itL}_{12}({\itx})-{\itr}({\itx})',-1)

В результате получаем, что погрешность интерполирования вблизи границ промежутка при |x|>0.72 растет с ростом степени интерполяционного полинома.

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

Рассмотрим еще  один пример, исследованный Бернштейном  в 1916г. Бернштейн рассматривал приближение  функции ƒ(x)=|x| на отрезке [-1,1] интерполяционными полиномами с равноотстоящими узлами и показал, что в этом случае не будет даже поточечной сходимости ни в одной точке отрезка [-1,1], кроме ±1.

Приблизим функцию  ƒ(x)=|x| на отрезке [-1,1] интерполяционными полиномами с равноотстоящими узлами, последовательно увеличивая число узлов интерполирования, и построим график зависимости погрешности интерполирования ƒ(x)-L_n(x) для x=0.9 от степени интерполяционного полинома, т.е. график зависимости ƒ(0.9)-L_n(0.9) от n.

% задание вектора  с числом узлов интерполирования

K=2:20;

% задание пустого  массива для последующей записи  в него ошибки

ERR=[];

% интерполяция  полиномами в цикле по узлам

for k=K

    % задание  k равноотстоящих узлов интерполирования

    x=linspace(-1,1,k);   

    % вычисление  в них значения |x|

    y=abs(x);

    % вычисление  значений интерполяционного полинома  в точке 0.9

    yy=lagrange(x,y,0.9);

    % вычисление  ошибки интерполяции в точке  0.9

    err=0.9-yy;

    % добавление  ошибки в массив

    ERR=[ERR err];

end

% создание графического  окна и осей

figure('Color','w')

axes

hold on

% вывод зависимости  ошибки от степени интерполяционного  полинома

plot(K-1,ERR)

В результате получаем, что в точке x=0.9 интерполяционный полином L_n(x) не сходится к функции ƒ(x)=|x|, т.е. нет даже поточечной сходимости, следовательно равномерной сходимости также нет.

Приведем несколько  классических теорем о сходимости интерполяционных полиномов, в которых имеется  ввиду, что интерполяция осуществляется на бесконечной последовательности, состоящей из наборов узлов интерполирования (интерполяционных сеток):

Здесь верхний  индекс обозначает номер интерполяционной сетки, а нижний - номер узла.

На первый взгляд кажется, что чем больше узлов  в сетке с увеличением ƒ(x), тем все ближе и ближе интерполяционный полином будет к интерполируемой функции. Однако, как утверждает теорема Фабера, это не так.

Теорема (Фабер). 
Для любой последовательности интерполяционных сеток найдется некоторая непрерывная функция ƒ(x), для которой последовательность соответствующих интерполяционных полиномов не сходится равномерно к ƒ(x).

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

Теорема (Марцинкевич). 
Для каждой непрерывной функции ƒ(x) существует последовательность интерполяционных сеток, такая, что построенные по этой последовательности интерполяционные полиномы равномерно сходятся к функции ƒ(x).

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

Теорема.

Если ƒ(x) целая функция, то тогда последовательность интерполяционных полиномов, построенных по любой последовательности сеток с узлами из отрезка [a,b], сходится равномерно к ƒ(x) на [a,b].

Следует заметить, что существование всех производных  у функции ƒ(x) недостаточно для того, чтобы утверждение этой теоремы было верным. Обычно приводят пример, выбирая в качестве интерполируемой функции следующую кусочно-заданную функцию на отрезке [-1,1]:

Очевидно, что  функция ƒ(x) непрерывна вместе со всеми своими производными на отрезке [-1,1] (да и вообще для любого x). Однако, если брать такую последовательность интерполяционных сеток, что все ее узлы принадлежат отрезку [-1,0], то тогда всегда L_n(x) ≡ 0 и для любого x > 0 никакой сходимости не будет.

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

Теорема. 
Если у функции ƒ(x) существует ограниченная на отрезке [a,b] производная ƒ'(x), а последовательность интерполяционных сеток состоит из наборов чебышевских узлов, то последовательность соответствующих интерполяционных полиномов равномерно сходится к ƒ(x) на отрезке [a,b].

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

1. В.В. Иванов. Методы вычислений на ЭВМ. Справочное пособие. Изд-во "Наукова думка". Киев. 1986.
2. Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. Численные методы. Изд-во "Лаборатория базовых знаний". 2003.
3. И.С. Березин, Н.П. Жидков. Методы вычислений. Изд. ФизМатЛит. Москва. 1962.
4. К. Де Бор. Практическое руководство по сплайнам. Изд-во "Радио и связь". Москва. 1985.
5. Дж. Форсайт, М.Мальком, К. Моулер. Машинные методы математических вычислений. Изд-во "Мир". Москва. 1980.