Будування інтерполяційних поліномів і їх обчислення

http://www..ru/.

МІНІСТЕРСТВО ОСВІТИ І НАУКИ, МОЛОДІ ТА СПОРТУ УКРАЇНИ МІЖНАРОДНИЙ НАУКОВО-ТЕХНІЧНИЙ УНІВЕРСИТЕТ ІМЕНІ АКАДЕМІКА ЮРІЯ БУГАЯ.

Контрольна робота.

з дисципліни «Чисельні методи».

Виконав:

Студент дистанційної форми навчання.

Групи ЗКд-31.

Кирилюк Богдан Олександрович.

КИЇВ — 2015.

Практична робота № 1.

X0=1.1 Y0=1.03×1=2.4.

X1=1.6 Y1=1.39×2=3.5.

X2=2.1 Y2=1.65×3=1.14.

X3=2.6 Y3=1.80.

X4=3.1 Y4=1.85.

X5=3.6 Y5=1.82.

Нехай наша функція f (х) задана на відрізку [a, b] =[1.1; 3.6] вузлами інтерполяції.

Необхідно побудувати інтерполяційні поліноми і обчислити з їх допомогою значення f (x_i) для заданих x_i.

x1=2.4×2=3.5×3=1.14.

1. Введемо два вектори і заповнимо їх початковими даними:

>> X = [1.1 1.6 2.1 2.6 3.1 3.6];.

>> y = [1.03 1.39 1.65 1.80 1.85 1.82];.

Нанесемо точки нашої табличної функції маркерами на графік за допомогою функції plot:

>> plot (X, y, 'ko').

Точкові значення табличної функції.

Рис. 1.

Тепер розглянемо роз’язування задачі інтерполяції цієї табличної функції різними методами..

1. Інтерполяція поліномами за методом найменших квадратів.

Стандартний метод інтерполяції, реалізований в MATLAB — метод найменших квадратів. Ми знаємо, що в MATLAB поліноми представляються у вигляді вектора своїх коефіцієнтів. Побудова поліномів фіксованого ступеня для наближення табличної функції однієї змінної робиться за допомогою функції polyfit з трьома вхідними параметрами.

polyfit=(х, у, n).

де х — вектор значень абсцис;

увектор значень ординат;

n-ступінь полінома.

Функція polyfit дозволяє знаходити коефіцієнти полінома.

Побудуємо інтерполяційний поліном 4-го ступеня для нашого прикладу і знайдемо його коефіцієнти:

>> P4 = polyfit (X, y, 4).

P4 = 0.0100 -0.0844 0.0472 0.9578 0.0169.

Для обчислення значення полінома в заданій точці (або точках) призначена відома нам функція polyval, у якості аргументів якої указуються вектор коефіцієнтів і вектор значень x_i.

Задамо діапазон значень x_i у вигляді вектора t (на відрізку від 1.0 до 3.6 з кроком 0.05:

t=1.0:0.05:3.6.

>> t=1.0:0.05:3.6.

Обчислимо вектор P4 значень полінома в цих точках.

>> p4=polyval (P4, t).

Побудуємо графік полінома за допомогою функції plot у вигляді суцільної лінії, причому виведемо його в те саме вікно, що і точковий графік. Додамо заголовок на графіці:

>> hold on.

>> plot (t, p4, 'k-').

>> hold on.

Наближення табличної функції поліномом 4-ступеня.

Рис. 2.

Знайдемо за допомогою побудованого інтерполяційного полінома значення g (x_i) для заданих x_i.

x1=2.4×2=3.5×3=1.14.

>> y1 = polyval (P4, 2.4).

>> y2 = polyval (P4, 3.5).

>> y3 = polyval (P4, 1.14).

Отримані результати представимо в таблиці:

Таблиця 1.

Обчислимо значення інтерполяційного полінома g (x_i) в заданих точках вектора.

X = [1.1 1.6 2.1 2.6 3.1 3.6];

>> P = polyval (P4,X).

P = 1.0300 1.3902 1.6496 1.8004 1.8498 1.8200.

Порівнюючи значення, отримані за допомогою інтерполяційного полінома 4-го ступеню із заданими в прикладі значеннями y_k, можна зробити висновок, що при заданій точності обчислень вони співпадають.

y = [1.03 1.39 1.65 1.80 1.85 1.82];

Підсумуємо виконані дії у вигляді файл-програми interpolation1.m.

%файл interpolation1.m.

%Функцію задано таблицею х-вектор вузлів інтерполяції.

% в — вектор значень функції у вузлах.

X = [1.1 1.6 2.1 2.6 3.1 3.6];

y = [1.03 1.39 1.65 1.80 1.85 1.82];

%будуємо точковий графік табличної функції.

plot (X, y, 'ko').

grid on.

%знаходимо коефіцієнти полінома 4-го ступеня.

P4=polyfit (X, y, 4);

%знаходимо точкові значення полінома на відрізку t.

t=1.0:0.05:3.6;

p4=polyval (P4,t);

%затримка.

hold on.

%будуємо графік полінома P4 по точках t в тому ж вікні.

plot (t, p4, 'k-').

hold on.

%оформляємо графік підписами.

title ('Approaching table function 4-degree polynomial.').

%Обчислимо значення поліному в заданих точках.

y1=polyval (P4, 2.4).

y2=polyval (P4, 3.5).

y3=polyval (P4, 1.14).

і виконаємо його Результати.

y1 = 1.7530.

y2 = 1.8308.

y3 = 1.0620.

2. Інтерполяційний поліном Лагранжа.

Інтерполяція методом Лагранжа припускає обчислення невідомих значень функції шляхом отримання середньозваженого значення функції у відомих сусідніх точках..

Поліном Лагранжа має наступний загальний вигляд:

або в іній формі запису:

L_n(х)=₀(х)f (х₀)+ ₁(х)f (х₁+.+ _k-1(х)f (x_k-1).

Коефіцієнти інтерполяційного полінома Лагранжа обчислюються за наступною формулою:

Таблиця 2.

• Необхідно побудувати інтерполяційний поліном Лагранжа і обчислити з його допомогою значення f (x_i) для заданих x_i.
• x1=2.4×2=3.5×3=1.14

1. Знайдемо коефіцієнти полінома w_i, і=0,5..

w0=((х-1.6)*(х-2.1)*(х-2.6)*(х-3.1)*(х-3.6))/((1.1−1.6)*(1.1−2.1)*(1.1−2.6)*(1.1−3.1)*(1.1−3.6)).

w1=((х-1.1)*(х-2.1)*(х-2.6)*(х-3.1)*(х-3.6))/((1.6−1.1)*(1.6−2.1)*(1.6−2.6)*(1.6−3.1)*(1.6−3.6)).

w2=((х-1.1)*(х-1.6)*(х-2.6)*(х-3.1)*(х-3.6))/((2.1−1.1)*(2.1−1.6)*(2.1−2.6)*(2.1−3.1)*(2.1−3.6)).

w3=((х-1.1)*(х-1.6)*(х-2.1)*(х-3.1)*(х-3.6))/((2.6−1.1)*(2.6−1.6)*(2.6−2.1)*(2.6−3.1)*(2.6−3.6)).

w4=((х-1.1)*(х-1.6)*(х-2.1)*(х-2.6)*(х-3.6))/((3.1−1.1)*(3.1−1.6)*(3.1−2.1)*(3.1−2.6)*(3.1−3.6)).

w5=((х-1.1)*(х-1.6)*(х-2.1)*(х-2.6)*(х-3.1))/((3.6−1.1)*(3.6−1.6)*(3.6−2.1)*(3.6−2.6)*(3.6−3.1)).

Всі рівняння коефіцієнтів містять змінну х, значення якої невідоме. Підставляючи в ці рівняння значення х1=2.4, х2=3.5 і х3=1.14, отримаємо чисельні значення коефіцієнтів в цих точках.

Призначимо змінній х значення першої точки х1=2.4.

>> x=2.4.

Введемо рівняння коефіцієнтів у вікні команд і отримаємо:

>> w0=((x-1.6)*(x-2.1)*(x-2.6)*(x-3.1)*(x-3.6))/((1.1−1.6)*(1.1−2.1)*(1.1−2.6)*(1.1−3.1)*(1.1−3.6)).

w0 = 0.0108.

>> w1=((x-1.1)*(x-2.1)*(x-2.6)*(x-3.1)*(x-3.6))/((1.6−1.1)*(1.6−2.1)*(1.6−2.6)*(1.6−3.1)*(1.6−3.6)).

w1 = -0.0874.

>> w2=((x-1.1)*(x-1.6)*(x-2.6)*(x-3.1)*(x-3.6))/((2.1−1.1)*(2.1−1.6)*(2.1−2.6)*(2.1−3.1)*(2.1−3.6)).

w2 = 0.4659.

>> w3=((x-1.1)*(x-1.6)*(x-2.1)*(x-3.1)*(x-3.6))/((2.6−1.1)*(2.6−1.6)*(2.6−2.1)*(2.6−3.1)*(2.6−3.6)).

w3 = 0.6989.

>> w4=((x-1.1)*(x-1.6)*(x-2.1)*(x-2.6)*(x-3.6))/((3.1−1.1)*(3.1−1.6)*(3.1−2.1)*(3.1−2.6)*(3.1−3.6)).

w4 = -0.0998.

>> w5=((x-1.1)*(x-1.6)*(x-2.1)*(x-2.6)*(x-3.1))/((3.6−1.1)*(3.6−1.6)*(3.6−2.1)*(3.6−2.6)*(3.6−3.1)).

w5 = 0.0116.

>> Тепер обчислимо значення полінома Лагранжа в цій точці.

P (х)=₀(х)f (х₀)+ ₁(х)f (х₁)+.+ _k-1(х)f (x_k-1).

де f (xi)=y_i.

Таблиця 3.

• >> P=w0*1.03+w1*1.39+w2*1.65+w3*1.80+w4*1.85+w5*1.82
• P = 1.7529

Таким же чином обчислюємо коефіцієнти і значення полінома в інших точках.

х2=3.5.

w0 = 0.0255.

w1 = -0.1613.

w2 = 0.4378.

w3 = -0.6810.

w4 = 0.7661.

w5 = 0.6129.

P = 1.8314.

Таким же чином обчислюємо коефіцієнти і значення полінома в інших точках.

х3=1.14.

w0 = 0.8290.

w1 = 0.3604.

w2 = -0.3454.

w3 = 0.2271.

w4 = -0.0846.

w5 = 0.0135.

P = 1.0618.

Результати подамо в таблиці.

Таблиця 4.

Отримані значення практично співпадають з отриманими за методом найменших квадратів !

Таблиця 5.

Узагальнимо метод Лагранжа для знаходження значень в будь-яких точках заданого відрізка. Напишемо файл-функцію Lagrange. m і тестову програму TestLagran.m.

Файл lagrange.m.

function yy=lagrange (x, y, xx).

% обчислення інтерполяційного полінома у формі Лагранжа.

% х — масив координат вузлів.

% в — масив значень інтерпольованої функції.

% xx — масив значень аргументу, для яких треба обчислити значення полінома.

% yy — масив значень полінома в точках xx.

% обчислюємо число вузлів інтерполяції (N=n+1).

N=length (x);

% створюємо нульовий масив значень інтерполяційного полінома.

yy=zeros (size (xx));

% в циклі рахуємо суму по вузлах.

for k=1:N.

% обчислюємо добутки.

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.

Файл TestLagran.m.

% TestLagran.m.

x = [1.1 1.6 2.1 2.6 3.1 3.6];

y = [1.03 1.39 1.65 1.80 1.85 1.82];

xx=[2.4 3.5 1.14];

[yy]= lagrange (x, y, xx).

% побудова графіків.

% виведення вузлів інтерполяції(х).

plot (x, y, 'ko').

hold on.

% виведення графіка полінома.

plot (xx, yy,'r').

Виконуємо файл TestLagran. m у Матлабі та отримуємо.

>> TestLagran.

yy = 1.7529 1.8314 1.0618.

Рис. 4.

3. Інтерполяційний поліном Ньютона.

Інтерполяційний поліном Ньютона ступеня _k, y_k) можна представити в наступній загальній формі:

Pn (х)=d_0,0+d_1,1(х-х₀)+…+d_{n, n}(х-х₀)(х-х₁)(х-x_n-1).

де d_{k, k}= d[х₀, х₁,…, x_k] - скінчені різниці.

d_{k, 0} = y_k і d_{k, j}= (d_{k, j-1} — d_k-1,j-1)/(x_k — x_k-j).

Для побудови і обчислення полінома Ньютона напишемо функцію [C, D] = newpoly (х, у) і збережемо її у файлі newpoly.m.

Файл newpoly.m.

function [C, D]=newpoly (x, y).

%Вхід х — вектор абсцис.

% y — вектор ординат.

%Вихід C — вектор коефіцієнтів полінома Ньютона.

% D — таблиця скінчених різниць.

n=length (x);

D=zeros (n, n);

D (, 1)=y'; %транспонування рядка в стовпець.

%Формування таблиці скінчених різниць D.

for j=2:n.

for k=j:n.

D (k, j)=(D (k, j-1)-D (k-1,j-1))/(x (k)-x (k-j+1));

end.

end.

% Обчислення коефіцієнтів полінома Ньютона.

C=D (n, n);

for k=(n-1):-1:1.

C=conv (C, poly (x (k)));

m=length (C);

C (m)=C (m)+D (k, k);

End.

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

Файл Testnewpoly.m.

%Тест TestNewpoly.

x=[1.1 1.6 2.1 2.6 3.1 3.6];

y=[1.03 1.39 1.65 1.80 1.85 1.82];

[C, D]=newpoly (x, y).

%Значения полінома в заданих точках.

y1=polyval (C, 2.4).

y2=polyval (C, 3.5).

y3=polyval (C, 1.14).

Виконаємо TestNewpoly в робочому середовищі Matlab.

>> TestNewpoly.

y1 = 1.7529.

y2 = 1.8314.

y3 = 1.0618.

4. Інтерполяція сплайнами.

Найпростішим способом інтерполяції є наближення даних сплайном нульового порядку (на кожній ділянці ступінь полінома дорівнює нулю), при якому значення в кожній проміжній точці приймається рівним найближчому значенню, заданому в таблиці. В результаті дані наближаються східчастою функцією, а саме наближення називається інтерполяцією по сусідніх точках. Лінійна інтерполяція заснована на з'єднанні сусідніх точок відрізками прямих — табличні дані наближаються ламаною лінією. Для отримання більш гладкої функції слід застосовувати інтерполяцію кубічними сплайнами. Всі ці способи реалізовані у функції Matlab interpl, у якості вхідних аргументів якої указуються координати абсцис і ординат табличної функції, координати абсцис проміжних точок, в яких обчислюються значення полінома і спосіб інтерполяції. Текстові константи, які задають спосіб інтерполяції:

`nearest' - наближення по сусідніх елементах;

`linear' - лінійна інтерполяція;

`spline' - інтерполяція кубічними сплайнами.

Вихідним аргументом interp1 є вектор значень полінома в проміжних точках.

Напишемо файл-програму, яка демонструє всі ці способи інтерполяції і обчислює значення функції в заданих точках.

Файл interplSpline.m.

x=[1.1 1.6 2.1 2.6 3.1 3.6];

y=[1.03 1.39 1.65 1.80 1.85 1.82];

%проміжні точки інтерполяції.

xi=[x (1):0.01:x (length (x))];

%задані точки інтерполяції.

xx=[2.4 3.5 1.14];

% побудова графіків.

% виведення вузлів інтерполяції(х) на графік.

plot (x, y,'ko').

hold on.

% обчислення східчастої функції в проміжних точках.

ynear=interp1(x, y, xi,'nearest');

% обчислення кусково-лінійної функції в проміжних точках.

yline=interp1(x, y, xi,'linear');

% обчислення кубічного сплайна в проміжних точках.

yspline=interp1(x, y, xi,'spline');

plot (xi, ynear,'k', xi, yline,'k:', xi, yspline,'k-.').

%title ('Способи інтерполяції функцій').

xlabel ('itx').

ylabel ('ity').

%legend ('таблична функція','по сусіднім елементам','лінійна','кубічними сплайнами').

%Обчислення значень функції в заданих точках.

ynear=interp1(x, y, xx,'nearest').

yline=interp1(x, y, xx,'linear').

yspline=interp1(x, y, xx,'spline').

Виконуємо файл interplSpline. m у Матлабі та отримуємо значення та малюнок.

>> interplSpline.

ynear = 1.8000 1.8200 1.0300.

yline = 1.7400 1.8260 1.0588.

yspline = 1.7529 1.8314 1.0621.

Рис. 7 'Способи інтерполяції функцій'.

Представимо знайдені результати в таблиці.

Таблиця 6.

поліном інтерполяція диференціювання функція.

Практична робота № 2.

Диференціювання полінома Нехай функція f (х) задана на відрізку [a, b] =[1,1.5] вузлами інтерполяції, поданими в таблиці 7.

Таблиця 7.

• Необхідно побудувати інтерполяційні поліноми і обчислити з їх допомогою значення f'(xi) в заданих точках xi.

х1=2.4.

х2=3.5.

х3=1.14.

1. Введемо два вектори і заповнимо їх початковими даними:.

х=[1.1 1.6 2.1 2.6 3.1 3.6];.

у=[1.03 1.39 1.65 1.80 1.85 1.82];.

Як ми знаємо, в MATLAB поліноми представляються у вигляді вектора своїх коефіцієнтів. Скористаємося функцією polyfit=(х, у, n) для знаходження коефіцієнтів полінома 4-го ступеню..

polyfit=(х, у, n).

де х — вектор значень абсцис;.

увектор значень ординат;.

n-ступінь полінома..

>> x = [1.1 1.6 2.1 2.6 3.1 3.6];.

>> y = [1.03 1.39 1.65 1.80 1.85 1.82];.

>> p4=polyfit (x, y, 4);.

Відповідь: >> p4.

p4 = 0.0100 -0.0844 0.0472 0.9578 0.0169.

Якщо n — ступінь полінома, то аналітичний запис його буде наступним:.

0.0100x⁴ — -0.0844x³ + 0.0472x² — 0.9578x + 0.0169.

Для обчислення похідної полінома призначена функція polyder. Виклик цієї функції з одним аргументом — вектором коефіцієнтів полінома, приводить до обчислення вектора коефіцієнтів похідної полінома..

Знайдемо коефіцієнти вектора першої похідної полінома р4.

>> d1 = polyder (p4).

Відповідь:.

d1 = 0.0400 -0.2531 0.0944 0.9578.

Отже, ми отримали коефіцієнти полінома 3-й ступеню:.

0.0400×3 — 0.2531×2 + 0.0944x + 0.9578.

Переконатися в правильності отриманої відповіді можна «вручну» продиференціювавши поліном р4..

Знайдемо значення похідної в заданих точках:.

>> df1 = polyval (d1, 2.4).

Відповідь:.

df1 = 0.2794.

>> df2 = polyval (d1, 3.5).

Відповідь:.

df2 = -0.0973.

>> df3 = polyval (d1, 1.14).

Відповідь:.

df3 = 0.7957.

Знайдемо похідні вищих порядків..

>> d2 = polyder (d1).

d2 = 0.1200 -0.5062 0.0944.

>> d3 = polyder (d2).

d3 = 0.2400 -0.5062.

>> d4 = polyder (d3).

d4 = 0.2400.

Геометричний зміст похідної.

Побудуємо графік інтерполяційного полінома на відрізку інтерполяції..

t=1.0:0.05:1.5.

>> t=1.0:0.05:3.6.

>> P4 = polyval (p4, t).

Рис. 8.

У тому самому вікні побудуємо графік дотичної в т. х1=1.26. Похідна.

>> df1 = polyval (d1,2.4).

df1 = 0.2794.

Рівняння дотичної, яка проходить через задану точку x0:.

y-y₀=y'(x-x₀).

Для нашого прикладу це буде рівняння:.

Y= 1.7530+0.2794*(t-2.4).

>> Y = 1.7530+0.2794*(t-2.4).

Рис. 9 Наближення табличної функції та дотична Скінчені різниці.

Диференціал функції дорівнює df (x) = f'(х)dx.

dx_i = xi = x_i+1-x_i — прирост аргументу.

df (x_i)= yi = y_i+1 — y_i — прирост функції.

Геометричний зміст диференціала Диференціал функції графічно зображається приростом ординати дотичної, проведеної до заданої точки..

Прирости функцій називаються скінчені різниці..

Знайдемо різниці між точками нашої табличної функції.

Таблиця 8.

• Побудуємо таблицю різниць її значень
• Таблиця 9

	x1= x2-x1.	x2=x3-x2.	x3=x4-x3.	x4=x5-x4.	x5=x6-x5.
xi.	0.1.	0.1.	0.1.	0.1.	0.1.
	y1= y2-y1.	y2=y3-y2.	y3=y4-y3.	y4=y5-y4.	y5=y6-y5.
yi.	0.02.	0.03.	0.03.	0.04.	0.05.

Якщо відстані між вузлами однакові x1=x2=x3 = h., то доцільно побудувати різниці другого та третього порядків для значень функції ²y_k, ³y_k..

Таблиця 10.

2y1= y2-y1.	2y2=y3-y2.	2y3=y4-y3.	2y4=y5-y4.
3y1=2y2 — 2y1.	3y2=2y3 — 2y3.	3y3=2y4 — 2y3.

Для знаходження вектора різниць вказаного порядку в Matlab призначена функція diff, яка може викликатися з одним або двома аргументами:.

dx=diff (x) — обчислення різниць 1-го порядку,.

dⁿx=diff (x, n) — обчислення різниць n-го порядку..

Знайдемо різниці 1-го та 2-го порядку для векторів табличної функції та диференціали dy/dx d²y/dx в цих точках..

Напишемо програму myDif. m у вигляді файл-програми..

% myDif. m — знаходження диференціала табличної функції.

Файл myDif. m.

x=[1.1 1.6 2.1 2.6 3.1 3.6];.

y=[1.03 1.39 1.65 1.80 1.85 1.82];.

dx=diff (x).

dy=diff (y).

df1=dy./dx.

d2y=diff (y, 2).

Отримаємо значення за допомогою нашої функції nyDif.

>> myDif.

Результати обчислень запишемо в таблицю.

Таблиця 11.

x.	dx.	y.	dy.	dy/dx.	d2y.
1.1.	0.5.	1.03.	0.36.	0.72.	— 0.1.
1.6.	0.5.	1.39.	0.26.	0.52.	— 0.11.
2.1.	0.5.	1.65.	0.15.	0.3.	— 0.1.
2.6.	0.5.	1.80.	0.05.	0.1.	— 0.08.
3.1.	0.5.	1.85.	— 0.03.	— 0.06.
3.6.		1.82.

Диференціювання полінома Ньютона Таблиці різниць використовуються для знаходження похідної табличної функції, для наближення якої використовується інтерполяційний поліном Ньютона..

Поліном Ньютона має наступний загальний вигляд:.

P_n(x)=d_0,0+d_1,1(x-x₀)+…+d_{n, n}(x-x₀)(x-x₁)…(x-x_n-1),.

где d_{k, k}= d[x₀, x₁,…x_k] - скінчені різниці..

d_{k, 0} = y_k и d_{k, j}= (d_{k, j-1} — d_k-1,j-1)/(x_k — x_k-j).

Для знаходження похідної функції в точках x_i (1.

f'(x_i) = (y_i-1 + y_i)/2)/h.

Ця формула непридатна для диференціювання в граничних точках відрізку. Для цього скористаємося наступними формулами другого порядку:.

f'(x₁) = (y₁ — ²y₁)/2)/h.

f'(x_n) = (y_n-1 — ²y_n-2)/2)/h.

Файл-програма myDiff. m демонструє використання різниць для знаходження першої похідної у вузлових точках..

Файл myDiff. m.

x=[1.1 1.6 2.1 2.6 3.1 3.6];.

y=[1.03 1.39 1.65 1.80 1.85 1.82];.

n=length (x);.

h=x (2)-x (1);.

dif=zeros (n);.

dx=diff (x);.

dy=diff (y);.

df1=dy./dx;.

d2y=diff (y, 2);.

dif (1)=((dy (1)-d2y (1))/2)/h;.

dif (6)=((dy (5)-d2y (4))/2)/h;.

for k=2:(n-1).

dif (k)=((dy (k-1)+ dy (k))/2)/h;.

end.

dif.

Запускаємо програму і отримаємо:.

>> myDiff.

Крім цього, для знаходження похідної інтерполяційного полінома Ньютона можна використати раніше розроблену функцію newpoly..

Скористаємося функцією [C, D] = newpoly (x, y), яка знаходить коефіцієнти інтерполяційного полінома Ньютона..

function [C, D]=newpoly (x, y).

%Вхід x — вектор абсцис.

% y — вектор ординат.

%Вихід C — вектор коефіцієнтів полінома Ньютона.

% D — таблиця різниць.

n=length (x);.

D=zeros (n, n);.

D (, 1)=y'; %транспонування рядка в столбець.

%Формування таблиці різниць D.

for j=2:n.

for k=j:n.

D (k, j)=(D (k, j-1)-D (k-1,j-1))/(x (k)-x (k-j+1));.

end.

end.

%Знаходження коефіцієнтів полінома Ньютона.

C=D (n, n);.

for k=(n-1):-1:1.

C=conv (C, poly (x (k)));.

m=length (C);.

C (m)=C (m)+D (k, k);.

end.

Тест-програма TestdiffNew, в якій задаються вхідні дані, викликається функція newpoly, обчислюються значення полінома і похідної в заданих точках..

Файл TestdiffNew. m.

%Тест TestdiffNew.

x=[1.1 1.6 2.1 2.6 3.1 3.6];.

y=[1.03 1.39 1.65 1.80 1.85 1.82];.

%Знаходження коефіцієнтів полінома Ньютона і таблиці скінечних різниць.

[C, D]=newpoly (x, y);.

%Знаходження коефіцієнтів першої похідної полінома Ньютона.

d1=polyder (C).

%Знаходження значення першої похідної у вузлових точках.

dif= polyval (d1,x).

%Знаходження значення першої похідної в заданих точках.

df1= polyval (d1,2.4).

df2=polyval (d1,3.5).

df3=polyval (d1,1.14).

Запускаємо файл-програму.

>> TestdiffNew.

Результат виконання програми :.

Вектор коефіцієнтів похідної полінома Ньютона:.

d1 = -0.0133 0.1653 -0.6788 0.7109 0.6382.

Значення похідної у вузлових точках:.

dif = 0.7993 0.6277 0.4093 0.1943 0.0127 -0.1257.

Значення похідної в заданих точках:.

df1 = 0.2776.

df2 = -0.1013.

df3 = 0.7888.

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

Практична робота № 3.

Мета роботи: Чисельне інтегрування функцій Знаходження невизначених інтегралів для табличних функцій. Обчислення визначених інтегралів для табличних функцій. Інтегрування функцій, заданих аналітично..

1 Знаходження невизначених інтегралів для табличних функцій.

Якщо підінтегральна функція задана поліномом, то для її інтегрування можна скористатися функцією.

q=polyint (p, c).

де c — константа інтегрування (вільний член первісної функції). Якщо с не задано, то вважається за 0.

>> q=polyint (p).

Функція повертає вектор коефіцієнтів первісної.

Приклад.

Нехай функція f (х) задана на відрізку [a, b] =[1,3.6] вузлами інтерполяції, поданими в таблиці.

Таблиця 12.

xi.	1.1.	1.6.	2.1.	2.6.	3.1.	3.6.
yi.	1.03.	1.39.	1.65.	1.80.	1.85.	1.82.

1. Знайдемо її інтерполяційний поліном.

Введемо два вектори і заповнимо їх початковими даними:

x=[1.1 1.6 2.1 2.6 3.1 3.6].

y=[1.03 1.39 1.65 1.80 1.85 1.82].

Знайдемо вектори полінома.

p4=polyfit (x, y,4).

>> p4=polyfit (x, y,4).

Відповідь:

p4 = 0.0100 -0.0844 0.0472 0.9578 0.0169.

Якщо n — ступінь полінома, то аналітичний запис його буде наступним:

0.0100x⁴ — 0.0844x³ + 0.0472x² + 0.9578x + 0.0169.

2. Знайдемо інтеграл полінома за допомогою функції polyint.

q=polyint (p4).

>> q=polyint (p4).

q = 0.0020 -0.0211 0.0157 0.4789 0.0169 0.

Аналітичний запис первісної.

q = 0.002x⁵ — 0.0211x⁴ + 0.0157x³ + 0.4789x² + 0.0169x.

Якщо тепер зайти похідну від q, то отримаємо наш поліном р4.

Функцію polyint можна використовувати для знаходження невизначеного інтеграла, тобто для відтворення первісної функції від табличної.

2. Обчислення визначених інтегралів для табличних функцій Візьмемо ту саму функцію f (х) яка задана на відрізку [a, b] =[1,3.6] вузлами інтерполяції.

Таблиця 13.

xi.	1.1.	1.6.	2.1.	2.6.	3.1.	3.6.
yi.	1.03.	1.39.	1.65.	1.80.	1.85.	1.82.

• Обчислимо інтеграл табличної функції за допомогою квадратурних формул та порівняємо результат.
• Введемо два вектори і заповнимо їх початковими даними:
• x=[1.1 1.6 2.1 2.6 3.1 3.6]
• y=[1.03 1.39 1.65 1.80 1.85 1.82]

Знайдемо відрізок h=0.1.

Обчислення за формулою прямокутників Файл my_rectangle.m.

function my_rectangle.

%h=(b-a)/(n-1).

x=[1.1 1.6 2.1 2.6 3.1 3.6];

y=[1.03 1.39 1.65 1.80 1.85 1.82];

n=6;

h=(x (n) — x (1))/(n-1);

sy=0;

for i=1:n-1.

sy=sy+y (i);

end.

In=h*sy.

Виконаємо програму в робочому вікні Matlab.

>> my_rectangle.

In = 3.8600.

Обчислення за формулою трапецій.

Для цього можна скористатися функцією trapz (x, y).

>> x=[1.1 1.6 2.1 2.6 3.1 3.6];

>> y=[1.03 1.39 1.65 1.80 1.85 1.82];

>> In = trapz (x, y).

Результат:

In = 4.0575.

Можна написати просту файл програму, яка обчислює формулу трапецій Файл my_trapz.m.

%my_trapz.m.

%h=(b-a)/n крок інтегрування.

x=[1.1 1.6 2.1 2.6 3.1 3.6];

y=[1.03 1.39 1.65 1.80 1.85 1.82];

n=6;

h=(x (n) — x (1))/(n-1);

sy=0;

for i=2:n-1.

sy=sy+y (i);

end.

In=h*((y (n)+y (1))/2+sy).

Результат виконання програми:

>> my_trapz.

In = 4.0575.

Знайдений розв’язок не відрізняється від отриманого за допомогою функції trapz (x, y).

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

Для табличної функції залишковий член формули трапецій можна обчислити за наближеною формулою:

Rтр = -((b-a)²/12)mean (²y).

де mean — середнє значення на відрізку інтегрування,.

²y — вектор скінчених різниць другого порядку.

В Matlab для обчислення скінчених різниць використовується функція diff з двома вхідними аргументами: вектором значень функції і порядку різниці.

Обчислимо залишковий член формули трапецій для нашого прикладу.

Виконаємо :

>> d2y=diff (y, 2).

Отримаємо :

d2y = -0.1000 -0.1100 -0.1000 -0.0800.

Виконаємо :

>> Rtr=-((0.5²)/12)* mean (d2y).

Отримаємо :

Rtr = 0.0020.

Додамо обчислене значення залишкового члена до обчисленого інтеграла:

>> In+Rtr.

Уточнене значення інтегралу:

ans = 4.0595.

3 Інтегрування функцій, заданих аналітично.

Інтерполяційний поліном нашої табличної функції має вигляд.

0.0100x⁴ — 0.0844x³ + 0.0472x² + 0.9578x + 0.0169.

syms x;

P=0.0100*x⁴ — 0.0844*x³ + 0.0472*x² + 0.9578*x + 0.0169.

Виконуємо у Матлабі :

>> P=0.01*x.^4 — 0.0844*x.^3 + 0.0472*x.^2 + 0.9578*x + 0.0169.

Отримуємо значення P.

P = 1.0299 1.3900 1.6493 1.7998 1.8488 1.8185.

>>Ip=int (P).

>> Ip=double (P);

>> Ip.

Ip = 1.0299 1.3900 1.6493 1.7998 1.8488 1.8185.

Висновок. Порівняємо значення інтеграла на відрізку [1.1,3.6], отримані різними способами Таблиця 13.

Спосіб.	Значення інтегралу.
За формулою прямокутників.	3.8600.
За формулою трапецій з урахуванням похибки обчислень.	4.0595.
З допомогою функції trapz (x, y).	4.0575.
З допомогою функції int (P, 1.1,3.6).	Не знайшов.

Практична робота № 4.

Розв’язування диференційних рівнянь Знайти розв’язок диференційного рівняння першого порядку [0,1].

y'=cos (1.5x+y)+(x-y), y (0)=0.

Розв’язок Оскільки це рівняння першого порядку, немає потреби у приведенні до системи рівнянь. Просто запишемо рівняння у вигляді допоміжної функції.

Файл oscil2. m:

function Dif=oscil2(x, y).

Dif=[cos (1.5*x+y)+(x-y)];

Для знаходження розв’язку рівняння скористаємося солвером ode113. Напишемо файл-функцію difur2.

файл difur2.m.

function difur2.

y0=[0];

[X, Y]=ode113(@oscil2,[0,1], y0).

plot (X, Y (, 1),'r.-').

hold on.

xlabel ('x').

ylabel ('y').

Виконаємо файл difur2 у командному рядку Matlab.

>> difur2.

X =.

0.0000.

0.0000.

0.0001.

0.0001.

0.0002.

0.0005.

0.0010.

0.0020.

0.0040.

0.0081.

0.0162.

0.0324.

0.0648.

0.1295.

0.2295.

0.3295.

0.4295.

0.5295.

0.6295.

0.7295.

0.8295.

0.9295.

1.0000.

Y =.

0.0000.

0.0000.

0.0001.

0.0001.

0.0002.

0.0005.

0.0010.

0.0020.

0.0040.

0.0081.

0.0162.

0.0323.

0.0645.

0.1274.

0.2181.

0.2978.

0.3645.

0.4175.

0.4575.

0.4861.

0.5050.

0.5162.

0.5205.

Результати запишемо в таблицю і виведемо на графік Таблиця наближених значень інтегралу рівняння.

y'=cos (1.5x+y)+(x-y), y (0)=0, x= [0,1].

Таблиця 14.

x.	y.
0.0000.	0.0000.
0.0000.	0.0000.
0.0001.	0.0001.
0.0001.	0.0001.
0.0002.	0.0002.
0.0005.	0.0005.
0.0010.	0.0010.
0.0020.	0.0020.
0.0040.	0.0040.
0.0081.	0.0081.
0.0162.	0.0162.
0.0324.	0.0323.
0.0648.	0.0645.
0.1295.	0.1274.
0.2295.	0.2181.
0.3295.	0.2978.
0.4295.	0.3645.
0.5295.	0.4175.
0.6295.	0.4575.
0.7295.	0.4861.
0.8295.	0.5050.
0.9295.	0.5162.
1.0000.	0.5205.

Рис. 10.