Титульный лист
Лабораторная работа по теме "Построение интерполяционных многочленов Эрмита "
Оглавление
Задание. 3
Постановка задачи. 3
Анализ задания. 3
Текст программы.. 4
скрипт lb5.m.. 4
скрипт Hpoly. m.. 5
Формулы найденных полиномов. 6
График функций. 6
Анализ результатов. 6
Задание
Построить таблицу разделенных разностей и интерполяционные многочлены Эрмита для функции
взяв в качестве узлов интерполяции x0=0 и x1=10.
Интерполяционные многочлены построить для следующих трех случаев:
1) Узел x0 двухкратный, узел x1 однократный;
2) Узел x0 однократный, узел x1 двухкратный;
3) Узел x0 двухкратный, узел x1 двухкратный;
В отчет включить постановку задачи, таблицы разделенных разностей, формулы полиномов, для каждого случая, графики исходной и интерполирующих функций. Для каждого из3 случаев оценить погрешность аппроксимации.
Постановка задачи
Согласно известным алгоритмам строим полиномы Эрмита, оцениваем погрешность и строим соответствующие графики.
Анализ задания
Функцию исходную задаем при помощи скрипт-файла f. m. График исходной функции строится с ее использованием – задается вектор x и вектор y абсцисс и ординат исходной функции.
В качестве вектора узлов интерполяции сначала задаем вектор ksi1 [0 0 10]. Для этого вектора находим ему соответствующий вектор y1, значения которого – значения функции f(x) в точках ksi1. Также задаем вектор значений первой производной в них, dy1. Затем пользуемся отдельной функцией Hpoly строим на базе значений функции и производной полином Эрмита третей степени. Сама функция написана так, чтобы работала не на интервале из двух точек, а на интервале, разбитом на большее число точек, для примера график такой интерполяции приведен в конце отчета. Здесь же в функции строится кривая, задаваемая полиномом, и выводятся значения коэффициентов полинома на график:
Погрешность оценивается как максимальное абсолютное отклонение найденного полинома от исходной функции.
Аналогичные действия повторяем и для других случаев.
Графики строим обычным способом, задавая несколько линий на одном объекте figure, а также оформляя график с использованием операндов legend, xlabel, ylabel, title.
Также в программе приведен пример построения многочлена Эрмита встроенной функцией Pchip.
Текст программы
скрипт lb5.m
% создание графического окна и осей
figure('Color','w')
axes
title('Интерполяция Эрмита при разной кратности узлов')
ksi1=[0 0 10];
y1 = f(ksi1);
disp('Интерполяция встроенной функцией pchip по узлам 0 5 10')
% вывод данных маркерами
subplot(2,2,1);
plot(ksi1,y1,'ro')
hold on;
x=0:0.1:10;
y=f(x);
plot(x, y,'-b');
hold on
yy = pchip([0 5 10], f([0 5 10]), x);
plot(x, yy,'--g');
disp('максимальная абсолютная погрешность')
e0=max(yy-y)
legend('точки интерполяции','исходная f(x)','интерполяция Эрмита pchip',4)
subplot(2,2,2);
disp('Интерполяция по узлам 0 0 10')
dy1=12.5*exp(-0.5*ksi1)./(5*exp(-0.5*ksi1)+1).^2;
plot(ksi1,y1,'ro')
hold on;
Hpoly(ksi1,y1,dy1)
hold on;
legend('исходная f(x)','x_0 ^(^2^) x_1 ^(^1^)', 4)
subplot(2,2,3);
disp('Интерполяция по узлам 0 10 10')
ksi2=[0 10 10];
y2 = f(ksi2);
plot(ksi2,y2,'ro')
hold on;
dy2=12.5*exp(-0.5*ksi2)./(5*exp(-0.5*ksi2)+1).^2;
Hpoly(ksi2,y2,dy2)
hold on;
legend('исходная f(x)','x_0 ^(^1^) x_1 ^(^2^) ', 4)
subplot(2,2,4);
disp('Интерполяция по узлам 0 0 10 10')
ksi3=[0 0 10 10];
y3 = f(ksi3);
plot(ksi3,y3,'ro')
hold on;
dy3=12.5*exp(-0.5*ksi3)./(5*exp(-0.5*ksi3)+1).^2;
Hpoly(ksi3,y3,dy3)
hold on;
legend('исходная f(x)','x_0 ^(^2^) x_1 ^(^2^)', 4)
скрипт Hpoly. m
function Hpoly(x, y,y1)
% определение числа узлов
n=length(x);
% нахождение всех коэффициентов a(k) и b(k)
a = y;
b = y1;
% вычисление разделенных разностей
DivDif = diff(y)./diff(x);
d=diff(x);
for i=1:length(DivDif)
if d(i)==0
DivDif(i)=y1(i)/1;
end
end
disp('вектор значений x');
x
disp('Таблица разделенных разностей')
DivDif
% вычисление всех коэффициентов c(k) и d(k)
for k = 1:n-1
c(k) = (3*DivDif(k) - 2*y1(k) - y1(k+1))/(x(k+1)-x(k));
d(k) = (y1(k) - 2*DivDif(k) + y1(k+1))/(x(k+1)-x(k))^2;
end
% построение графиков каждого из полиномов
for k=1:n-1
% вычисление k-го полинома в 30 точках между x(k) и x(k+1)
xx=linspace(x(k),x(k+1),30);
yy=a(k) + b(k)*(xx-x(k)) + c(k)*(xx-x(k)).^2 + d(k)*(xx-x(k)).^3;
err(k)=max(abs(yy-f(xx)));
% выбор случайного цвета
rgbcolor = [rand rand rand];
% построение графика k-го полинома
plot(xx, yy,'Color', rgbcolor,'Linestyle','-')
% вывод коэффициентов k-го полинома
text('Position',[xx(15) yy(15)],...
'String',num2str([a(k); b(k); c(k); d(k)]), 'Color', rgbcolor,...
'VerticalAlignment','middle','HorizontalAlignment','center')
end
disp('максимальная абсолютная погрешность')
max(err)
Формулы найденных полиномов
Ниже приведены найденные программой таблицы разделенных разностей и погрешности интерполяционных полиномов для каждого случая.
Интерполяция встроенной функцией pchip по узлам 0 5 10
максимальная абсолютная погрешность
e0 =
0.3525
Интерполяция по узлам 0 0 10
вектор значений x
x =
0 0 10
Таблица разделенных разностей
DivDif =
0.3472 0.4004
максимальная абсолютная погрешность
ans =
0.3769
Интерполяция по узлам 0 10 10
вектор значений x
x =
0 10 10
Таблица разделенных разностей
DivDif =
0.4004 0.0788
максимальная абсолютная погрешность
ans =
0.3769
Интерполяция по узлам 0 0 10 10
вектор значений x
x =
0 0 10 10
Таблица разделенных разностей
DivDif =
0.3472 0.4004 0.0788
максимальная абсолютная погрешность
ans =
0.3769
Коэффициенты полиномов приведены на графике.
Таким образом, сведем найденную информацию в таблицу
Разделенная разность | Полином | |||
узлы | 1 случай | 2случай | 3 случай |
|
1-2 | 0.3472 | 0.4004 | 0.3472 | |
2-3 | 0.4004 | 0.0788 | 0.4004 | |
3-4 | 0.0788 | |||
погрешность | 0.3769 | 0.3769 | 0.3769 |
График функций
Графики исходной заданной функции, и построенных полиномов для разных случаев. Впрочем из графиков видно, что полиномы совпадают.
Пример работы программы если диапазон от 0 до 10 разбить с шагом 2. Каждые две точки соединены интерполяционным полиномом Эрмита третей степени.
Анализ результатов
Согласно примеру можно построить полиномы Эрмитапо несколько другой схеме.
Результат второго расчета:
Интерполяция встроенной функцией pchip по узлам 0 5 10
максимальная абсолютная погрешность
e0 =
0.3525
Интерполяция по узлам 0 0 10
вектор значений ksi1
ksi1 =
0 0 10
Таблица разделенных разностей
DivDif =
0.3472 0.4004
DivDif2 =
0.0053
максимальная абсолютная погрешность
Коэффициенты полинома: 0.83333
Коэффициенты полинома: 0.34722
Коэффициенты полинома:0.0053149
Интерполяция по узлам 0 10 10
вектор значений ksi2
ksi2 =
0 10 10
Таблица разделенных разностей
DivDif =
0.4004 0.0788
DivDif2 =
-0.0322
максимальная абсолютная погрешность
Коэффициенты полинома: 0.83333
Коэффициенты полинома: 0.40037
Коэффициенты полинома:-0.032155
Интерполяция по узлам 0 0 10 10
вектор значений ksi3
ksi3 =
0 0 10 10
Таблица разделенных разностей
DivDif =
0.3472 0.4004 0.0788
DivDif2 =
0.0053 -0.0322
DivDif3 =
-0.0037
максимальная абсолютная погрешность
Коэффициенты полинома: 0.83333
Коэффициенты полинома: 0.34722
Коэффициенты полинома:0.0053149
Коэффициенты полинома:-0.003747
По этим данным можно построить таблицы разделенных разностей, полиномы и указать погрешности
узлы | Разности | Полином | Погрешность | ||
1 порядка | 2 порядка | 3 порядка | |||
1 случай | |||||
1-2 | 0.3472 | 0.0053 |
| 0.8958 | |
2-3 | 0.4004 | ||||
2 случай | |||||
1-2 | 0.4004 | -0.0322 |
| 3.2155 | |
2-3 | 0.0788 | ||||
2 случай | |||||
1-2 | 0.3472 | 0.0053 |
| 3.7470 | |
2-3 | 0.4004 | -0.0322 | |||
3-4 | 0.0788 |
Графики показаны на рисунке ниже. Видно, что такой подход имеет меньшую точность, чем предложенный выше. Вторая схема представлена в скрипте lb52.m.
Скрипт
% создание графического окна и осей
figure('Color','w')
axes
title('Интерполяция Эрмита при разной кратности узлов')
ksi1=[0 0 10];
y1 = f(ksi1);
disp('Интерполяция встроенной функцией pchip по узлам 0 5 10')
% вывод данных маркерами
subplot(2,2,1);
plot(ksi1,y1,'ro')
hold on;
x=0:0.1:10;
y=f(x);
plot(x, y,'-b');
hold on
yy = pchip([0 5 10], f([0 5 10]), x);
plot(x, yy,'--g');
disp('максимальная абсолютная погрешность')
e0=max(yy-y)
legend('точки интерполяции','исходная f(x)','интерполяция Эрмита pchip',4)
subplot(2,2,2);
disp('Интерполяция по узлам 0 0 10')
dy1=12.5*exp(-0.5*ksi1)./(5*exp(-0.5*ksi1)+1).^2;
plot(ksi1,y1,'ro')
hold on;
%построение полинома
% определение числа узлов
n=length(ksi1);
% нахождение всех коэффициентов
a = y1(1);
% вычисление разделенных разностей
DivDif = diff(y1)./diff(ksi1);
d=diff(ksi1);
for i=1:length(DivDif)
if d(i)==0
DivDif(i)=dy1(i)/1;
end
end
disp('вектор значений ksi1');
ksi1
disp('Таблица разделенных разностей')
DivDif
DivDif2=diff(DivDif)./diff([0 10])
b=DivDif(1);
c=DivDif2(1);
xx=linspace(ksi1(1),ksi1(n),30);
yy=a + b*(xx-ksi1(1)) + c*(xx-ksi1(1)).^2;
disp('максимальная абсолютная погрешность')
err=max(abs(yy-f(xx)))
% выбор случайного цвета
rgbcolor = [rand rand rand];
% построение графика k-го полинома
plot(xx, yy,'Color', rgbcolor,'Linestyle','-')
% вывод коэффициентов k-го полинома
disp(strcat('Коэффициенты полинома:',num2str([a; b; c])))
%полином построен
hold on;
legend('исходная f(x)','x_0 ^(^2^) x_1 ^(^1^)', 4)
subplot(2,2,3);
disp('Интерполяция по узлам 0 10 10')
ksi2=[0 10 10];
y2 = f(ksi2);
plot(ksi2,y2,'ro')
hold on;
dy2=12.5*exp(-0.5*ksi2)./(5*exp(-0.5*ksi2)+1).^2;
%построение полинома
% определение числа узлов
n=length(ksi2);
% нахождение всех коэффициентов
a = y2(1);
% вычисление разделенных разностей
DivDif = diff(y2)./diff(ksi2);
d=diff(ksi2);
for i=1:length(DivDif)
if d(i)==0
DivDif(i)=dy2(i)/1;
end
end
disp('вектор значений ksi2');
ksi2
disp('Таблица разделенных разностей')
DivDif
DivDif2=diff(DivDif)./diff([0 10])
b=DivDif(1);
c=DivDif2(1);
xx=linspace(ksi2(1),ksi2(n),30);
yy=a + b*(xx-ksi2(1)) + c*(xx-ksi2(1)).^2;
disp('максимальная абсолютная погрешность')
err=max(abs(yy-f(xx)))
% выбор случайного цвета
rgbcolor = [rand rand rand];
% построение графика k-го полинома
plot(xx, yy,'Color', rgbcolor,'Linestyle','-')
% вывод коэффициентов k-го полинома
disp(strcat('Коэффициенты полинома:',num2str([a; b; c])))
%полином построен
hold on;
legend('исходная f(x)','x_0 ^(^1^) x_1 ^(^2^) ', 4)
subplot(2,2,4);
disp('Интерполяция по узлам 0 0 10 10')
ksi3=[0 0 10 10];
y3 = f(ksi3);
plot(ksi3,y3,'ro')
hold on;
dy3=12.5*exp(-0.5*ksi3)./(5*exp(-0.5*ksi3)+1).^2;
%построение полинома
% определение числа узлов
n=length(ksi3);
% нахождение всех коэффициентов
a = y3(1);
% вычисление разделенных разностей
DivDif = diff(y3)./diff(ksi3);
d=diff(ksi3);
for i=1:length(DivDif)
if d(i)==0
DivDif(i)=dy3(i)/1;
end
end
disp('вектор значений ksi3');
ksi3
disp('Таблица разделенных разностей')
DivDif
DivDif2=diff(DivDif)./diff([0 10])
DivDif3=diff(DivDif2)./diff([0 10])
b=DivDif(1);
c=DivDif2(1);
dc=DivDif3(1);
xx=linspace(ksi3(1),ksi3(n),30);
yy=a + b*(xx-ksi3(1)) + c*(xx-ksi3(1)).^2+dc*(xx-ksi3(1)).^3;
disp('максимальная абсолютная погрешность')
err=max(abs(yy-f(xx)))
% выбор случайного цвета
rgbcolor = [rand rand rand];
% построение графика k-го полинома
plot(xx, yy,'Color', rgbcolor,'Linestyle','-')
% вывод коэффициентов k-го полинома
disp(strcat('Коэффициенты полинома:',num2str([a; b; c; dc])))
%полином построен
hold on;
legend('исходная f(x)','x_0 ^(^2^) x_1 ^(^2^)', 0)
