Будування інтерполяційних поліномів і їх обчислення
Реалізація інтерполяції поліномами за методами найменших квадратів і Лагранжа в Matlab. Наближення даних сплайном нульового порядку. Диференціювання полінома. Геометричний зміст похідної. Чисельне інтегрування функцій. Розв’язування диференційних рівнянь.
| Рубрика | Программирование, компьютеры и кибернетика |
| Вид | контрольная работа |
| Язык | украинский |
| Дата добавления | 01.06.2015 |
| Размер файла | 285,3 K |
Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже
Студенты, аспиранты, молодые ученые, использующие базу знаний в своей учебе и работе, будут вам очень благодарны.
Размещено на http://www.allbest.ru/
МІНІСТЕРСТВО ОСВІТИ І НАУКИ, МОЛОДІ ТА СПОРТУ УКРАЇНИ
МІЖНАРОДНИЙ НАУКОВО-ТЕХНІЧНИЙ УНІВЕРСИТЕТ ІМЕНІ АКАДЕМІКА ЮРІЯ БУГАЯ
Контрольна робота
з дисципліни “Чисельні методи”
Виконав:
Студент дистанційної форми навчання
Групи ЗКд-31
Кирилюк Богдан Олександрович
КИЇВ - 2015
Практична робота №1
X0=1.1 Y0=1.03 x1=2.4
X1=1.6 Y1=1.39 x2=3.5
X2=2.1 Y2=1.65 x3=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(xi) для заданих xi.
x1=2.4 x2=3.5 x3=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, у якості аргументів якої указуються вектор коефіцієнтів і вектор значень xi.
Задамо діапазон значень xi у вигляді вектора 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(xi) для заданих xi.
x1=2.4 x2=3.5 x3=1.14
>> y1 = polyval(P4, 2.4)
>> y2 = polyval(P4, 3.5)
>> y3 = polyval(P4, 1.14)
Отримані результати представимо в таблиці:
Таблиця 1
|
xi |
2.4 |
3.5 |
1.14 |
|
|
yi |
1.7530 |
1.8308 |
1.0620 |
Обчислимо значення інтерполяційного полінома g(xi) в заданих точках вектора
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-го ступеню із заданими в прикладі значеннями yk, можна зробити висновок, що при заданій точності обчислень вони співпадають.
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. Інтерполяційний поліном Лагранжа
Інтерполяція методом Лагранжа припускає обчислення невідомих значень функції шляхом отримання середньозваженого значення функції у відомих сусідніх точках.
Поліном Лагранжа має наступний загальний вигляд:
або в іній формі запису:
Ln(х)=0(х)f(х0)+ 1(х)f(х1+.+ k-1(х)f(xk-1)
Коефіцієнти інтерполяційного полінома Лагранжа обчислюються за наступною формулою:
Таблиця 2
|
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 |
- Необхідно побудувати інтерполяційний поліном Лагранжа і обчислити з його допомогою значення f(xi) для заданих xi.
- x1=2.4 x2=3.5 x3=1.14
1. Знайдемо коефіцієнти полінома wi, і=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(х)=0(х)f(х0)+ 1(х)f(х1)+.+ k-1(х)f(xk-1)
де f(xi)=yi
Таблиця 3
|
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 |
- >> 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
|
xi |
1.24 |
3.5 |
1.14 |
|
|
yi |
1.7529 |
1.8314 |
1.0618 |
Отримані значення практично співпадають з отриманими за методом найменших квадратів !
Таблиця 5
|
xi |
1.24 |
3.5 |
1.14 |
|
|
yi |
1.7530 |
1.8308 |
1.0620 |
Узагальнимо метод Лагранжа для знаходження значень в будь-яких точках заданого відрізка. Напишемо файл-функцію 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. Інтерполяційний поліном Ньютона
Інтерполяційний поліном Ньютона ступеня <n, який проходить через точки (xk, yk ) можна представити в наступній загальній формі:
Pn(х)=d0,0+d1,1(х-х0)+….+dn,n(х-х0)(х-х1)(х-xn-1)
де dk,k= d[х0,х1,....,xk] - скінчені різниці
dk,0 = yk і dk,j= (dk,j-1 - dk-1,j-1)/(xk - xk-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
|
Спосіб згладжування |
xi |
2.4 |
3.5 |
1.14 |
|
|
По найближчих точках (nearest) |
yi |
1.8000 |
1.8200 |
1.0300 |
|
|
Лінійна інтерполяція (linear) |
yi |
1.7400 |
1.8260 |
1.0588 |
|
|
Кубічним сплайном (spline) |
yi |
1.7529 |
1.8314 |
1.0621 |
поліном інтерполяція диференціювання функція
Практична робота №2
Диференціювання полінома
Нехай функція f(х) задана на відрізку [a,b] =[1,1.5] вузлами інтерполяції, поданими в таблиці 7
Таблиця 7
|
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 |
- Необхідно побудувати інтерполяційні поліноми і обчислити з їх допомогою значення 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.0100x4 - -0.0844x3 + 0.0472x2 - 0.9578x + 0.0169
Для обчислення похідної полінома призначена функція polyder. Виклик цієї функції з одним аргументом - вектором коефіцієнтів полінома, приводить до обчислення вектора коефіцієнтів похідної полінома.
Знайдемо коефіцієнти вектора першої похідної полінома р4
>> d1 = polyder(p4)
Відповідь:
d1 = 0.0400 -0.2531 0.0944 0.9578
Отже, ми отримали коефіцієнти полінома 3-й ступеню:
0.0400x3 - 0.2531x2 + 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-y0=y'(x-x0)
Для нашого прикладу це буде рівняння:
Y= 1.7530+0.2794*(t-2.4)
>> Y = 1.7530+0.2794*(t-2.4)
Рис. 9 Наближення табличної функції та дотична
Скінчені різниці
Диференціал функції дорівнює df(x) = f'(х)dx
dxi = xi = xi+1-xi - прирост аргументу
df(xi)= yi = yi+1 - yi - прирост функції
Геометричний зміст диференціала
Диференціал функції графічно зображається приростом ординати дотичної, проведеної до заданої точки.
Прирости функцій називаються скінчені різниці.
Знайдемо різниці між точками нашої табличної функції
Таблиця 8
|
x1 |
x2 |
x3 |
x4 |
x5 |
x6 |
||
|
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 |
- Побудуємо таблицю різниць її значень
- Таблиця 9
- Обчислимо інтеграл табличної функції за допомогою квадратурних формул та порівняємо результат.
- Введемо два вектори і заповнимо їх початковими даними:
- 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]
|
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.., то доцільно побудувати різниці другого та третього порядків для значень функції 2yk, 3yk.
Таблиця 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-го порядку,
dnx=diff(x,n) - обчислення різниць n-го порядку.
Знайдемо різниці 1-го та 2-го порядку для векторів табличної функції та диференціали dy/dx d2y/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 |
Диференціювання полінома Ньютона
Таблиці різниць використовуються для знаходження похідної табличної функції, для наближення якої використовується інтерполяційний поліном Ньютона.
Поліном Ньютона має наступний загальний вигляд:
Pn(x)=d0,0+d1,1(x-x0)+…+dn,n(x-x0)(x-x1)…(x-xn-1),
где dk,k= d[x0,x1,…xk] - скінчені різниці.
dk,0 = yk и dk,j= (dk,j-1 - dk-1,j-1)/(xk - xk-j)
Для знаходження похідної функції в точках xi (1<i<n) використовується формула (отримана шляхом диференціювання півсум формул Ньютона 1-го порядку):
f'(xi) = (yi-1 + yi)/2)/h
Ця формула непридатна для диференціювання в граничних точках відрізку. Для цього скористаємося наступними формулами другого порядку:
f'(x1) = (y1 - 2y1)/2)/h
f'(xn) = (yn-1 - 2yn-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.0100x4 - 0.0844x3 + 0.0472x2 + 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.002x5 - 0.0211x4 + 0.0157x3 + 0.4789x2 + 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 |
Знайдемо відрізок 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)2/12)mean(2y)
де mean - середнє значення на відрізку інтегрування,
2y - вектор скінчених різниць другого порядку.
В Matlab для обчислення скінчених різниць використовується функція diff з двома вхідними аргументами: вектором значень функції і порядку різниці.
Обчислимо залишковий член формули трапецій для нашого прикладу.
Виконаємо :
>> d2y=diff(y,2)
Отримаємо :
d2y = -0.1000 -0.1100 -0.1000 -0.0800
Виконаємо :
>> Rtr=-((0.5^2)/12)* mean(d2y)
Отримаємо :
Rtr = 0.0020
Додамо обчислене значення залишкового члена до обчисленого інтеграла:
>> In+Rtr
Уточнене значення інтегралу:
ans = 4.0595
3 Інтегрування функцій, заданих аналітично
Інтерполяційний поліном нашої табличної функції має вигляд
0.0100x4 - 0.0844x3 + 0.0472x2 + 0.9578x + 0.0169
syms x;
P=0.0100*x^4 - 0.0844*x^3 + 0.0472*x^2 + 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
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
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 |
0 |
|
|
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
Размещено на Allbest.ru
Подобные документы
Чисельне інтегрування, формула Сімпсона, значення інтегралу від функцій та формули трапецій. Знаходження коренів рівняння методом Ньютона. Наближення функцій поліномами вищого порядку. Метод Ейлера та його модифікації. Визначення похибок розрахунків.
контрольная работа [6,1 M], добавлен 04.07.2010Стандартний спосіб розв’язання задачі Коші для звичайного диференціального рівняння першого порядку чисельними однокроковими методами. Геометричний зміст методу Ейлера. Побудова графіку інтегральної кривої. Особливість оцінки похибки за методом Рунге.
курсовая работа [112,9 K], добавлен 30.11.2009Виконання "ручного" розв'язування рівняння методом Ньоютона. Розробка програми на мові С#, яка реалізує введення вихідних даних, розв'язання заданого рівняння, виведення результатів у зручній формі на екран. Визначення початкового наближення кореня.
лабораторная работа [120,9 K], добавлен 19.01.2022Вибір емпіричної формули. Метод оберненої матриці. Розв’язування систем лінійних рівнянь на ПК. Вибір двох апроксимуючих функцій. Розрахунки у середовищі MS Excel для лінійної функції, для квадратичної функції та у середовищі MS Visual Studio (мовою С#).
курсовая работа [658,8 K], добавлен 18.08.2014Загальні відомості та геометричний зміст розв'язання задачі Коші. Використання методу Ейлера для розв'язання звичайних диференціальних рівнянь першого порядку. Розробка блок-схеми та реалізація алгоритму в середовищі програмування Borland Delphi 7.0.
курсовая работа [398,1 K], добавлен 14.10.2012Структурна схема моделі (пакет MATLAB) та її описання. Математична модель у вигляді передавальних функцій, у вигляді диференційного рівняння. Алгоритм рішення (рекурентне співвідношення) та його програмна реалізація. Системи диференційних рівнянь.
курсовая работа [551,8 K], добавлен 14.02.2009Ортогонaлізування функцій. Порівняння дискретного та хвильового перетворення. Інтерполяційні поліноми Лагранжа і Ньютона. Метод найменших квадратів. Побудова кривої для заданих результатів вимірювань. Розв’язання задачі по Лапласу операційним методом.
курсовая работа [2,2 M], добавлен 10.04.2012Опис методів обчислення формули Ньютона-Котеса та поліномів Лежандра. Розгляд програмування процедур вводу меж інтегрування, ініціації елементів квадратурних формул Гауса та Чебишева. обчислення визначеного інтеграла і виводу результатів на екран.
курсовая работа [82,1 K], добавлен 23.04.2010Зародження системи Matlab. Високоефективна мова інженерних і наукових обчислень. Інтерактивна система, основним об'єктом якої є масив. Обчислення мінімумів, нулів функцій. Апроксимація й інтерполяція даних. Обчислення кінцевих різниць, перетворення Фур'є.
лабораторная работа [146,4 K], добавлен 18.01.2013Обґрунтування переваги чисельного диференціювання функції з використанням інтерполяційної формули Стірлінга по відношенню до формул Ньютона, Гауса та Бесселя. Розробка оптимального алгоритму обчислення другої похідної. Лістинг, опис і тестування програми.
курсовая работа [483,2 K], добавлен 21.10.2013
