">
Информатика Программирование
Информация о работе

Тема: Системы компьютерных вычислений

Описание: Вычисление интерполяционного полинома в форме Лагранжа. Ошибка интерполяции, чебышевские узлы. Сходимость интерполяционного полинома Ньютона . Нулевой массив и построение графика значений интерполяционного полинома. Результат работы программы на графике.
Предмет: Информатика.
Дисциплина: Программирование.
Тип: Курсовая работа
Дата: 08.08.2012 г.
Язык: Русский
Скачиваний: 14
Поднять уникальность

Похожие работы:

Московский авиационный институт

(национальный исследовательский университет)

Курсовая работа

Системы компьютерных вычислений

Вариант №….

Москва 2011/12 уч. Год

№1

С помощью функций A = randi(100*7+50,[10,10]) и f=randi(100*7+50,[10,1]) получим матрицы A и f.

A =

518 461 330 671 705 32 604 401 341 583

391 676 423 373 1 570 378 355 423 514

58 190 743 157 330 244 326 243 174 240

305 166 718 542 9 392 532 324 109 625

323 177 253 525 176 714 417 574 508 496

495 608 117 346 531 502 192 748 266 441

164 527 313 128 521 355 57 199 383 15

61 659 416 596 199 488 442 595 522 628

351 534 23 667 245 500 506 101 447 86

547 496 651 107 70 173 442 27 380 631

f =

254

323

596

107

457

698

478

36

524

193

u U U-u  n  0,3408 0,3411 0,0003  <0,001  0,1450 0,1451 0,0001  <0,001  0,1239 0,1240 0,0001  <0,001  1,6263 1,6278 0,0016  <0,01  0,4042 0,4046 0,0004  <0,001  0,7539 0,7546 0,0007  <0,001  2,3147 2,3168 0,0021  <0,01  0,9246 0,9254 0,0008  <0,001  0,4753 0,4757 0,0004  <0,001  1,3264 1,3276 0,0012  <0,01  

2) Решаем СЛАУ Au = f :

Для этого воспользуемся левым делением в MATLAB

>> u=Af соответствует A*u=f

3) Найдем абсолютную и относительную погрешности:

Решаем в MATLAB

>> A1=A.*0.0001

>> f1=f.*0.001

>> A2=A-A1

>> f2=f-f1

>> u1=A2f2

>> u2=abs(u-u1)

>>

U3=u2./u

^u=0.18% =0.076%

4) Округлим сомнительные цифры числа, оставив верные знаки:

(для удобства значения взяты по модулю):

№2

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

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

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

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

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

При программировании интерполяционного полинома по второму способу в функции 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;

Проинтерполируем данную нам функцию  на отрезке [-1,1.5] с шагом 0.5 и построим графики  и полученного интерполяционного полинома. В данном случае будет 7 узлов и, следовательно, интерполяционный полином получится 6-ой степени. Для построения графика интерполяционного полинома вычислим его значения при помощи функции lagrange в 1000 равномерно отстоящих друг от друга точках на отрезке [-1,1.5].

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

x=[-1;-1;-0.5;0;0.5;1;1.5];

y=((exp(-x))./sqrt(exp(-x)+1))-x;

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

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

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

yy=lagrange(x,y,xx);

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

figure(Color,w)

% вывод графика ((exp(-x))./sqrt(exp(-x)+1))-x

plot(y)

hold on

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

plot(xx,yy,r)

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

plot(x,y,bo)

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

legend(yitx,{itL_n}({itx}),nodes,-1)

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

x=[-1;-1;-0.5;0;0.5;1;1.5];

y=((exp(-x))./sqrt(exp(-x)+1))-x;

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

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

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

yy=lagrange(x,y,xx);

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

figure(Color,w)

% вывод графика ((exp(-x))./sqrt(exp(-x)+1))-x

fplot(((exp(-x))./sqrt(exp(-x)+1))-x, [-1;-1;-0.5;0;0.5;1;1.5])

hold on

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

plot(xx,yy,r)

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

plot(x,y,bo)

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

legend(yitx,{itL_n}({itx}),nodes,-1

Интерполяция функции ((exp(-x))./sqrt(exp(-x)+1))-x)полиномом 6-ой степени

Приближение полиномиальной функции. Выберем в качестве интерполируемой функции полином шестой степени:

P5(x)= -0.0004х5+0.0005х4-0.0170х3+0.1545х2-1.5302х+0.7071

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

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

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

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

p5=inline(-0.0004*x.^5+0.0005*x.^4-0.0170*x.^3+0.1545*x.^2-1.5302*x+0.7071);

% вычисление полинома 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)

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

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,p5(xx)-yy)

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

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

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



Интерполяция полинома шестой степени по шести узлам

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

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

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

P5(x)= -0.0004х5+0.0005х4-0.0170х3+0.1545х2-1.5302х+0.7071

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

p5=inline(-0.0004*x.^5+0.0005*x.^4-0.0170*x.^3+0.1545*x.^2-1.5302*x+0.7071)% задание трех равноотстоящих узлов на [-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)

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



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

Проинтерполируем полином

P5(x)= -0.0004х5+0.0005х4-0.0170х3+0.1545х2-1.5302х+0.7071

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

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

p5=inline(-0.0004*x.^5+0.0005*x.^4-0.0170*x.^3+0.1545*x.^2-1.5302*x+0.7071);

% задание десяти равноотстоящих узлов на [-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-ой степени), хотя она должна была бы равняться нулю (при вычислениях в точной арифметике). Заметим, что характер распределения ошибки сохраняется - чем ближе к границам отрезка интерполирования, тем ошибка больше.



Поведение ошибки для интерполяционных полиномов L9(x), L19(x), L29(x)

Вычислительные свойства интерполяционного полинома делают его неприменимым при достаточно большом числе узлов интерполирования. Изучим зависимость максимального уклонения интерполяционного полинома от интерполируемой функции снова на примере полинома P5(x)= -0.0004х5+0.0005х4-0.0170х3+0.1545х2-1.5302х+0.7071

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

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

p5=inline(-0.0004*x.^5+0.0005*x.^4-0.0170*x.^3+0.1545*x.^2-1.5302*x+0.7071);

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

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})

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



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

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

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

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



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

?n+1(x) = (x-x1)(x-x2)...(x-xn+1)

где x1,x2,...,xn+1 - узлы интерполяционного полинома Ln(x).

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

Продемонстрируем использование этой оценки для проверки ошибки интерполяции на примере интерполяции функции  на отрезке [-1,1.5] интерполяционным полиномом 5-ой степени с равномерно отстоящими семью узлами интерполяции. Для вычисления интерполяционного полинома L5(x) воспользуемся функцией lagrange:

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

x =

-1.0000 -0.5000 0 0.5000 1.0000 1.5000

>> y=((exp(-x))./sqrt(exp(-x)+1))-x

y =

2.4097 1.5130 0.7071 -0.0215 -0.6855 -1.2982

>> xx=pi/6

xx =

0.5236

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

yy =

-0.0542

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

err =

0.5542

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

est =

1.2450e-05

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

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

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

x=linspace(-11.5,7);

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

y=(((exp(-x))./sqrt(exp(-x)+1))-x);

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

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

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

yy=lagrange(x,y,xx);

err=abs(yy-(((exp(-xx))./sqrt(exp(-xx)+1))-(xx)));

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

est=zeros(size(xx));

for i=1:length(xx)

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

end

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

figure(Color,w)

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

legend(error,estimate)



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

?n(x) = (x-x1)(x-x2)...(x-xn+1) ? 0в узлах интерполяции x1,x2,...,xn+1.

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

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

?n+1(x) = (x-x1)(x-x2)...(x-xn+1) ? 0

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

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

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

Tk(x) = cos(k arccos x), |x|?1.

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

T1(x) = cos(arccos x) = x

В случае k = 2 тоже очевидно, что T2(x) есть полином второй степени, если воспользоваться известным тригонометрическим тождествомcos2? = 2cos?? - 1,положив ? = arccos x.

Тогда получаемT2(x) = 2x? - 1,

Тригонометрическое тождествоcos(k + 1)? = 2cos?cosk? - cos(k - 1)позволяет легко установить, что для полиномов Чебышева справедливо следующее рекуррентное соотношение :Tk+1(x) = 2xTk(x) - Tk-1(x)

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

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

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

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)



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

>> k=6/10/13;

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

>> y=(((exp(-x))./sqrt(exp(-x)+1))-x);

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

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

>> err=yy-sin(xx)

>> m=0;k-1;

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

>> ycheb=(((exp(-xcheb))./sqrt(exp(-xcheb)+1))-(xcheb));

>> yy=lagrange(xcheb,ycheb,xx);

>> errcheb=yy-sin(xx);

>> plot(xx,err,xx,errcheb)

>> legend(equidistance,chebyshev)

>> title(Interpolation error)



Ошибки интерполяции для 6 узлов



Ошибки интерполяции для 10 узлов



Ошибки интерполяции для 13 узлов

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

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

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

a=-1;

b=1.5;

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

fun=inline((((exp(-x))./sqrt(exp(-x)+1))-x));

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

k=6/10/13;

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

x=linspace(a,b,k);

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

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);

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

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)



Ошибки интерполяции для 6 узлов



Ошибки интерполяции для 10 узлов



Ошибки интерполяции для 13 узлов

Полиномы Ньютона















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



где



Мы можем считать, что точки пронумерованы в обратном порядке, тогда



где



Тогда при вычислении значения интерполяционного полинома в форме Ньютона требуется выполнить следующие действия:

найти разделенные разности, стоящие внизу треугольной таблицы

вычислить значение приведенного выше полинома от заданного значения x

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

Алгоритм вычисления разделенных разностей.



Цикл по 

Цикл по 



После завершения цикла



При выполнении второго пункта, т.е. при вычислении интерполяционного полинома по формуле



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



сложений и



умножений.

Записав выражение для интерполяционного полинома в виде



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

Для вычисления интерполяционного полинома в форме Ньютона в MATLAB можно воспользоваться приведенной ниже функцией newton

function yy = newton(x, y, xx)

%Вычисление интерполяционного полинома в форме Ньютона

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

% y – массив ординат точек, через которые должен проходить интерполяционный полином

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

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

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

N = length(x);

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

DIFF = y;

for k = 1 : N-1

for i = 1: N - k

DIFF(i) = (DIFF(i+1) - DIFF(i)) / (x(i+k) - x(i));

end

end

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

yy = DIFF(1) * ones(size(xx));

for k = 2 : N

yy = DIFF(k) + (xx - x(k)) .* yy;

end

Интерполируем при помощи функции newton следующую табличную функцию, заданную массивами

>> x=[-1 -0.5 0 0.5 1 1.5];

>> y=[2.4097 1.5130 0.7071 -0.0215 -0.6855 -1.2982];

при помощи интерполяционного полинома (5-ой степени) и вычислим значения интерполяционного полинома в 100 равноотстоящих точках на отрезке интерполирования [-1, 1.5]

xx = linspace(x(1), x(end), 100);

yy = newton(x, y, xx);

Для проверки работы функции newton можно построить график, отобразив значения табличной функции маркерами, а значения интерполяционного полинома непрерывной линией

plot(x,y,o,xx,yy)

В результате получаем

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

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

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



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



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

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

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



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

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

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

r=inline((((exp(-x))./sqrt(exp(-x)+1))-x));

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

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)


Функция Рунге и интерполяционные полиномы (сверху), ошибка интерполяции (снизу).

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

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

Приблизим функцию ?(x)=|x| на отрезке [-1,1] интерполяционными полиномами с равноотстоящими узлами, последовательно увеличивая число узлов интерполирования, и построим график зависимости погрешности интерполирования ?(x)-Ln(x) для x=0.9 от степени интерполяционного полинома, т.е. график зависимости ?(0.9)-Ln(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 интерполяционный полином Ln(x) не сходится к функции ?(x)=|x|, т.е. нет даже поточечной сходимости, следовательно равномерной сходимости также нет.


Зависимость ?(0.9)-Ln(0.9) от n для ?(x)=|x|.

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

В.В. Иванов. Методы вычислений на ЭВМ. Справочное пособие. Изд-во "Наукова думка". Киев. 1986.

Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. Численные методы. Изд-во "Лаборатория базовых знаний". 2003.

И.С. Березин, Н.П. Жидков. Методы вычислений. Изд. ФизМатЛит. Москва. 1962.

К. Де Бор. Практическое руководство по сплайнам. Изд-во "Радио и связь". Москва. 1985.

Дж. Форсайт, М.Мальком, К. Моулер. Машинные методы математических вычислений. Изд-во "Мир". Москва. 1980.

Джон Г. Мэтьюз, Куртис Д. Финк . Численные методы. Использование MATLAB. 2001.

Б.П. Демидович, И.А. Марон, Э.З. Шувалова. Численные методы анализа.1967.

Материалы с сайта www.matlab.exponenta.ru

Интернет-ресурсы:

http://эссе.рф - сборник не проиндексированных рефератов. Поиск по рубрикам и теме. Большинство текстов бесплатные. Магазин готовых работ.