МIНIСТЕРСТВО ОСВIТИ І НАУКИ УКРАЇНИ
Національний унiверситет "Львiвська полiтехнiка"
Кафедра САПР
Розрахунково–графічна робота
ЛЬВІВ 2014
1. МЕТА РОБОТИ
Мета роботи – закріпити знання отримані при вивченні курсу «Чисельні методи в інформатиці» та ознайомитися з новими відомими чисельними методами.
Варіант завдань:
Завдання 1
Розв'язати систему лінійних алгебраїчних рівнянь методом Гаусса за схемою Халецького наведену у Додатку 1 згідно варіанту індивідуального завдання отриманого у викладача.
Завдання 2
Для заданого варіанту обчислити значення функції в 5 різних точках для двох різних інтервалів інтерполяції. Дослідити вплив степені полінома (n1=1,2,3 n2=5,6,7 n3=8,9,10) на точність одержаних результатів.
Завдання 3
Хід роботиЗавдання 1
Текст програми на мові Matlab
clear
clc
a = [
-11 1 -6 -3 ;
-6 15 -4 3 ;
1 1 11 1;
-4 1 1 -11 ;
]
d = [ -9 ;4; -2; 4 ]
b = [
0 0 0 0;
0 0 0 0;
0 0 0 0;
0 0 0 0;
];
c = b;
n = 4;
for i = 1:n
b(i,i) = 1
b(1,i) = a(1,i)/a(1,1)
c(i,1) = a(i,1)
end
y(1) = d(1,1) / a(1,1)
for i = 2:n
sumc = 0;
sumb = 0;
sumy = 0;
for k = i:n
for j=1:i-1
sumc = sumc + c(k,j)*b(j,i)
end
c(k,i) = a(k,i) - sumc
end
for k = i+1:n
for j=1:i-1
sumb = sumb + c(i,j)*b(j,k);
end
b(i,k) = (a(i,k) - sumb)/c(i,i)
end
for j =1:i-1
sumy = sumy + c(i,j)*y(j);
end
y(i) = (d(i,1) - sumy)/c(i,i)
end
x(n) = y(n)
for i = 1:3
sumx = 0;
for j = n-i+1:n
sumx = sumx + b(n-i,j)*x(j);
end
x(n-i) = y(n-i) - sumx
end
Результати виконання програми
a =
-11 1 -6 -3
-6 15 -4 3
1 1 11 1
-4 1 1 -11
d =
-9
4
-2
4
b =
1 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
b =
1 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
c =
-11 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
b =
1 0 0 0
0 1 0 0
0 0 0 0
0 0 0 0
b =
1.0000 -0.0909 0 0
0 1.0000 0 0
0 0 0 0
0 0 0 0
c =
-11 0 0 0
-6 0 0 0
0 0 0 0
0 0 0 0
b =
1.0000 -0.0909 0 0
0 1.0000 0 0
0 0 1.0000 0
0 0 0 0
b =
1.0000 -0.0909 0.5455 0
0 1.0000 0 0
0 0 1.0000 0
0 0 0 0
c =
-11 0 0 0
-6 0 0 0
1 0 0 0
0 0 0 0
b =
1.0000 -0.0909 0.5455 0
0 1.0000 0 0
0 0 1.0000 0
0 0 0 1.0000
b =
1.0000 -0.0909 0.5455 0.2727
0 1.0000 0 0
0 0 1.0000 0
0 0 0 1.0000
c =
-11 0 0 0
-6 0 0 0
1 0 0 0
-4 0 0 0
y = 0.8182
sumc = 0.5455
c =
-11.0000 0 0 0
-6.0000 14.4545 0 0
1.0000 0 0 0
-4.0000 0 0 0
sumc = 0.4545
c =
-11.0000 0 0 0
-6.0000 14.4545 0 0
1.0000 0.5455 0 0
-4.0000 0 0 0
sumc = 0.8182
c =
-11.0000 0 0 0
-6.0000 14.4545 0 0
1.0000 0.5455 0 0
-4.0000 0.1818 0 0
b =
1.0000 -0.0909 0.5455 0.2727
0 1.0000 -0.0503 0
0 0 1.0000 0
0 0 0 1.0000
b =
1.0000 -0.0909 0.5455 0.2727
0 1.0000 -0.0503 0.5472
0 0 1.0000 0
0 0 0 1.0000
y = 0.8182 0.6164
sumc = 0.5455
sumc = 0.5180
c =
-11.0000 0 0 0
-6.0000 14.4545 0 0
1.0000 0.5455 10.4820 0
-4.0000 0.1818 0 0
sumc = -1.6638
sumc = -1.6730
c =
-11.0000 0 0 0
-6.0000 14.4545 0 0
1.0000 0.5455 10.4820 0
-4.0000 0.1818 2.6730 0
b =
1.0000 -0.0909 0.5455 0.2727
0 1.0000 -0.0503 0.5472
0 0 1.0000 0.0409
0 0 0 1.0000
y = 0.8182 0.6164 -0.3009
sumc = -1.0909
sumc = -0.9914
sumc = -0.8821
c =
-11.0000 0 0 0
-6.0000 14.4545 0 0
1.0000 0.5455 10.4820 0
-4.0000 0.1818 2.6730 -10.1179
y = 0.8182 0.6164 -0.3009 -0.7872
x = 0 0 0 -0.7872
x = 0 0 -0.2687 -0.7872
x = 0 1.0336 -0.2687 -0.7872
x = 1.2734 1.0336 -0.2687 -0.7872
Рішення СЛАР методом простої ітерації
Перш ніж застосовувати метод ітерацій, необхідно переставити рядки вихідної системи таким чином, щоб на діагоналі стояли найбільші за модулем коефіцієнти матриці. Якщо при цьому умова все таки не виконується, то іноді вдається забезпечити збіжність методу за допомогою наступного методу.
Нехай дана система Ax = b. Перетворимо її до вигляду: x = Qx + c
де Q = E - D • A, c = D • b
Тут D - деяка матриця. Нам необхідно підібрати таку матрицю D, щоб виконувалася умова | Q | <1.
Щоб отримати | Q | <1, використовуємо наступний спосіб.
Маємо СЛАР
A x = b (1)
Припускаючи, що aii ≠ 0 дозволимо нове рівняння системи (1) щодо x1, друга - щодо x2, ..., n-е рівняння - щодо xn. В результаті отримаємо:
x1 = β1 - α12x2 - α13x3 - ... - α1nxn
x2 = β2 - α21x1 - α23x3 - ... - α2nxn
xn = βn - αn1xn - αn3x3 - ... - αnn-1xn-1
де βi = bi / aii; αij = aij / aii при i ≠ j; αii = 0
Система (2) в матричній формі має вигляд:
x = β - αx
Систему будемо вирішувати методом послідовних наближень. Нехай x0 = β, тоді:
x1 = b - a x0
x2 = b - a x1
....
xk +1 = b - a xk
Перш ніж застосовувати метод, необхідно переставити рядки вихідної системи таким чином, щоб на діагоналі стояли найбільші за модулем коефіцієнти матриці.
-11
1
-6
-3
-6
15
-4
3
1
1
11
1
-4
1
1
-11
Наведемо до вигляду:
x1=0.82-0.0909x2+0.55x3+0.27x4x2=0.27-0.4x1-0.27x3+0.2x4x3=-0.18+0.0909x1+0.0909x2+0.0909x4x4=-0.36+0.36x1-0.0909x2-0.0909x3
Покажемо обчислення на прикладі кількох ітерацій.
N=1x1=0.82 - 0 • (-0.0909) - 0 • 0.55 - 0 • 0.27=0.82x2=0.27 - 0 • (-0.4) - 0 • (-0.27) - 0 • 0.2=0.27x3=-0.18 - 0 • 0.0909 - 0 • 0.0909 - 0 • 0.0909=-0.18x4=-0.36 - 0 • 0.36 - 0 • (-0.0909) - 0 • (-0.0909)=-0.36N=2x1=0.82 - 0.27 • (-0.0909) - (-0.18) • 0.55 - (-0.36) • 0.27=1.04x2=0.27 - 0.82 • (-0.4) - (-0.18) • (-0.27) - (-0.36) • 0.2=0.62x3=-0.18 - 0.82 • 0.0909 - 0.27 • 0.0909 - (-0.36) • 0.0909=-0.25x4=-0.36 - 0.82 • 0.36 - 0.27 • (-0.0909) - (-0.18) • (-0.0909)=-0.65N=3x1=0.82 - 0.62 • (-0.0909) - (-0.25) • 0.55 - (-0.65) • 0.27=1.19x2=0.27 - 1.04 • (-0.4) - (-0.25) • (-0.27) - (-0.65) • 0.2=0.75x3=-0.18 - 1.04 • 0.0909 - 0.62 • 0.0909 - (-0.65) • 0.0909=-0.27x4=-0.36 - 1.04 • 0.36 - 0.62 • (-0.0909) - (-0.25) • (-0.0909)=-0.71
Завдання 2.
Текст програми на мові Matlab
x21=[7 8 8.5 9 9.5 10];
x31=[7 7.5 8 8.5 9 9.5 9.75 10];
y21=exp(x21).*sin(x21);
y31=exp(x31).*sin(x31);
x=[6 7.75 8.75 9.25 11];
yt=exp(x).*sin(x);
x=[6 7.75 8.75 9.25 11];
z=(10);
for k=1:5
cs=0;
for j=1:6
p=1;
for i=1:6
if i~=j
p=((x(k)-x21(i))/(x21(j)-x21(i)))*p;
end
end
cs=cs+p*y21(j);
end
z(k)=cs;
end
x=[6 7.75 8.75 9.25 11];
h1=(10);
for k=1:5
cs=0;
for j=1:7
p=1;
for i=1:7
if i~=j
p=((x(k)-x31(i))/(x31(j)-x31(i)))*p;
end
end
cs=cs+p*y31(j);
end
h1(k)=cs;
end
abp11=abs(yt-y);
abp21=abs(yt-z);
abp31=abs(yt-h1);
vidp11=abp11./abs(yt);
vidp21=abp21./abs(yt);
vidp31=abp31./abs(yt);
figure (1)
hold on
grid on
plot(x,abp21,'.-')
plot(x,abp31,'--')
legend('pohybka pry n=5','pohybka pry n=6')
hold off
figure (2)
hold on
grid on
plot(x,vidp21,'.-')
plot(x,vidp31,'--')
legend('vidnosna pohybka pry n=5','vidnosna pohybka pry n=6')
hold off
figure (3)
plot(x,yt,'.-')
hold on
grid on
plot(x,z,'k:')
plot(x,h1,'--')
legend('Tochna funkcia','Nablyjhena funkcia pry n=5','Nablyjhena funkcia pry n=6')
hold off
x22=[10 11 12 13 14 15];
x32=[10 10.5 11 12.5 13 14 15];
y22=exp(x22).*sin(x22);
y32=exp(x32).*sin(x32);
x2=[9 10.75 11.5 13.5 16];
yt2=exp(x2).*sin(x2);
x2=[9 10.75 11.5 13.5 16];
z2=(10);
for k=1:5
cs=0;
for j=1:6
p=1;
for i=1:6
if i~=j
p=((x2(k)-x22(i))/(x22(j)-x22(i)))*p;
end
end
cs=cs+p*y22(j);
end
z2(k)=cs;
end
x2=[9 10.75 11.5 13.5 16];
h2=(10);
for k=1:5
cs=0;
for j=1:7
p=1;
for i=1:7
if i~=j
p=((x2(k)-x32(i))/(x32(j)-x32(i)))*p;
end
end
cs=cs+p*y32(j);
end
h2(k)=cs;
end
abp12=abs(yt2-y2);
abp22=abs(yt2-z2);
abp32=abs(yt2-h2);
vidp12=abp12./abs(yt2);
vidp22=abp22./abs(yt2);
vidp32=abp32./abs(yt2);
figure (4)
hold on
grid on
plot(x2,abp22,'.-')
plot(x2,abp32,'--')
legend('pohybka pry n=5','pohybka pry n=6')
hold off
figure (5)
hold on
grid on
plot(x,vidp22,'.-')
plot(x,vidp32,'--')
legend('vidnosna pohybka pry n=5','vidnosna pohybka pry n=6')
hold off
figure (6)
plot(x2,yt2,'.-')
hold on
grid on
plot(x2,z2,'k:')
plot(x2,h2,'--')
legend('Tochna funkcia','Nablyjhena funkcia pry n=5','Nablyjhena funkcia pry n=6')
hold off
Результати виконання програми
Для першого проміжку [7;10]
Абсолютні похибки:
Відносні похибки:
Графік точної та наближених функцій:
Для 2-го проміжку [10;15]
Абсолютні похибки:
Відносні похибки:
Графік точної та наближених функцій:
Завдання 3
Інтерполяційний поліном Ерміта.
У вузлах може бути задано значення функції f та похідні - f', f'', …, f(Ni-1) - усього Ni штук.
Тоді всього буде N0 + N1 + … + Nn умов.
Шукається поліном степеня
у якому Ni - кратність вузла і.
Поліном існує та єдиний. Кількість рівнянь дорівнює кількості коефіцієнтів. Розгляньмо поліном степеня m, похідні якого Hm(j)(xi) = 0, j = 0, 1, …, Ni . Це означає. що xi - корінь кратності Ni .
Тоді на [a; b] існує N0 + N1 + … + Nn = m+1 коренів. але степінь його m, тому при всіх x Hm(x) = 0.
Будемо шукати Hm(x) у вигляді
Hm(x) = Ln(x) + wn(x) Hm-n(x),
де
Ln(xi) = yi, wn(xi) = 0, i = 0, 1, …, n.
Похідна:
H'm(x) = L'n(x) + w'n(x) Hm-n(x) + wn(x) H'm-n(x),
причому
H'm(xi) = L'n(xi) + w'n(x) Hm-n(xi).
У тих точках. де задано f'(xi), знайдемо
Далі шукають другу похідну
H''m(x) = L''n(x) + w''n(x) Hm-n(x) + 2 w'n(x) H'm-n(x) + wn(x) H''m-n(x),
причому
H''m(xi) = L''n(xi) + w''n(xi) Hm-n(xi) + 2 w'n(xi) H'm-n(xi),
Звідки виражають H'm-n(xi), тому що інші величини є відомими, і т. ін..
Приклад.
Нехай задано значення функції та її похідних у точках x0 = 0; x1 = 1; x2 = 2:
f(0) = 1; f '(0) = 1; N0 = 2;
f(1) = 1,5; N1 = 1;
f(2) = 0; f '(2) = 0; f ''(2) = 1; N0 = 2;
Таким чином, m = (2 + 1 + 3) - 1 = 5. Поліном Ерміта шукаємо у вигляді H5(x) = L2(x) + w3(x) H2(x):
,
тому
w3'(x) = (x - 1)(x - 2) + x(x - 2) + x(x - 1);
w3'(0) = 2; w3'(1) = -1; w3'(2) = 2;
H5''(x) = -2 + w3''(x) H2(x) + 2w3'(x) H2'(x) + w3(x) H2''(x);
H5''(2) = -2 + w3''(2) H2(2) + 2w3'(2) H2'(2) + w3(2) H2''(2);
1 = -2 + 6(2-1) 5/4 + 4 H2''(2) + 0;
H2'(2) = -1,125. Далi є арифм. помилки...
H0(x) = const.
Програмна реалізація
#include <stdio.h>
#include <math.h>
#include <windows.h>
#include <iostream>
using namespace std;
double X,x[5],y[5],yp[5],Pn[5];
int i,n;
void main()
{
setlocale(LC_ALL,"Russian");
for(i=0;i<4;i++)
{
cout<<"Введите значение x["<<i<<"]: ";
cin>>x[i];
cout<<"Введите значение y["<<i<<"]: ";
cin>>y[i];
cout<<"Введите значение y'["<<i<<"]: ";
cin>>yp[i];
cout<<"-----------------------------------------"<<endl;
}
cout<<"Введите значение х для которого вычисляем значение функции: ";
cin>>X;
cout<<"-----------------------------------------"<<endl;
for(n=1;n<4;n++)
{
Pn[n]=y[0]+(X-x[0])*(yp[0]+(X-x[0])*(yp[0]-(y[0]-y[n])/(x[0]-x[n])+(X-x[n])*(yp[0]-2*(y[0]-y[n])/(x[0]-x[n])+y[n])/(x[0]-x[n]))/(x[0]-x[n]));
}
cout<<"Расчитывая по формуле:\n\nPn[n]=y[0]+(X-x[0])*(yp[0]+(X-x[0])*(yp[0]-(y[0]-y[n])/(x[0]-x[n])+\n+(X-x[n])*(yp[0]-2*(y[0]-y[n])/(x[0]-x[n])+y[n])/(x[0]-x[n]))/(x[0]-x[n]))\n\nзначения будут равны:"<<endl;
cout<<"P(n) на участке [x0,x1] ="<<Pn[1]<<endl;
cout<<"P(n) на участке [x0,x2] ="<<Pn[2]<<endl;
cout<<"P(n) на участке [x0,x3] ="<<Pn[3]<<endl;
}
Зведення задачі інтегрування системи звичайних диференціальних рівнянь до послідовності задач інтегрування диференціальних рівнянь першого порядку.
Інтегровні типи диференціальних рівнянь першого порядку, розв(язаних відносно похідної
Неповні рівняння.
а) Диференціальне рівняння, яке не містить шуканої функції
має вигляд
, . (2.33)
Припустимо, що f(x) являється неперервною функцією на . Тоді функція
(2.34)
є загальним розв`язком диференціального рівняння (2.33) в області
a < x < b, -< y < + . (2.35)
Особливих розв’язків диференціальне рівняння (2.33) не має.
Разом з диференціальним рівнянням (2.33) розглянемо початкові умови
. (2.36)
Проінтегруємо диференціальне рівняння (2.34) від до x
.
Знаходимо с з умови (2.36)
(2.37)
– загальний розв(язок диференціального рівняння (2.33) в формі Коші.
Якщо f(x) – неперервна на за виключенням точки , в якій приймає нескінченне значення, то замість диференціального рівняння (2.33) будемо розглядати рівняння
. (2.33()
Пряма є розв(язком диференціального рівняння (2.33() і ми цей розв(язок повинні приєднати до розв(язку диференціального рівняння (2.33). Цей розв(язок може бути частинним або особливим в залежності від того зберігається чи порушується в будь-якій його точці єдиність. Якщо – частинний розв(язок, то його часто можна отримати з загального при нескінченних значеннях с, якщо ж він є особливим, то його отримують з загального при .
Рівняння, яке не містить незалежної змінної має вигляд
. (2.38)
Припускаємо, що функція визначена і неперервна на інтервалі . Замість (2.38) розглянемо диференціальне рівняння
. (2.39)
Диференціальне рівняння (2.39) не містить шуканої функції і воно розв(язується аналогічно диференціальному рівнянню (2.33).
Якщо , y є (c,d), то
(2.40)
– загальний розв(язок диференціального рівняння (2.39) в області
c < y < d, -< x < + .
Аналогічно
(2.41)
– загальний інтеграл в формі Коші.
Якщо неперервна на (c,d) і приймає нульове значення при , то ми повинні розглядати диференціальне рівняння (2.38). Розв(язок буде частинним, якщо в кожній його точці зберігається єдиність і особливим, якщо в кожній його точці порушується єдиність. Якщо частинний розв(язок, то ми його отримуємо при нескінченних значеннях , якщо особливий, то при .
Якщо в точці перетворюється в нескінченність , то розглядаємо диференціальне рівняння (2.39), яке має неперервну праву частину на (c,d). При цьому диференціальне рівняння на має єдиний розв(язок (мал. 2.5).
Висновки
На даній графічно-розрахунковій роботі я розв’язав індивідуальну систему лінійних рівнянь методом Гауса–Халецького. Обчислив значення функції в 5 різних точках для двох різних інтервалах інтерполяції (Лагранж), дослідив вплив степенів полінома (n1=1,2,3 n2=5,6,7) на точність одержаних результатів, навів відповіді на теоретичні запитання.