Алгебраїчні методи розв'язання рівнянь
Обчислювальні методи розв’язку нелінійних рівнянь. Методи лінійної алгебри. Знаходження визначника матриці методом алгебраїчних доповнень. Інтерполювання функцій. Методи чисельного інтегрування функцій. Розв’язування звичайних диференціальних рівнянь.
| Рубрика | Математика |
| Вид | лекция |
| Язык | русский |
| Дата добавления | 13.09.2010 |
| Размер файла | 107,6 K |
Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже
Студенты, аспиранты, молодые ученые, использующие базу знаний в своей учебе и работе, будут вам очень благодарны.
68
ЗМІСТ
Обчислювальні методи розв'язку нелінійних рівнянь
Метод хорд
Метод дотичних
Метод хорд-дотичних
Метод ітерацій
Метод половинного поділу
МЕТОДИ ЛІНІЙНОЇ АЛГЕБРИ
Знаходження визначника матриці методом алгебраїчних доповнень
Розв'язок системи рівнянь методом Крамера
Розв'язок системи рівнянь методом Гаусса
Знаходження визначника матриці методом Гаусса
Розв'язок системи лінійних рівнянь за допомогою методу Зейделя
ІНТЕРПОЛЮВАННЯ ФУНКЦІЙ
Інтерполяційний многочлен Лагранжа
Інтерполяційний многочлен Ньютона
ЧИСЕЛЬНЕ ІНТЕГРУВАННЯ ФУНКЦІЙ
Метод прямокутників
Метод трапецій
Метод Сімпсона
РОЗВ'ЯЗУВАННЯ ЗВИЧАЙНИХ ДИФЕРЕНЦІАЛЬНИХ РІВНЯНЬ
Метод Ейлера
Уточнений метод Ейлера
Метод Рунге-Кутта 4-го порядку
Обчислювальні методи розв'язку нелінійних рівнянь
Нехай задано рівняння з однією змінною
f(x)=0, (1)
де функція f(x) визначена і неперервна на деякому проміжку <a;b>.
Розв'язати рівняння - означає знайти множину його коренів, тобто таких значень x<a;b>, при яких рівняння (1) перетворюється в тотожність. Корінь рівняння (1) називають ще нулем функції f(x). Якщо функція f(x) - алгебраїчний многочлен, то рівняння (1) називають алгебраїчним. Якщо функція f(x) містить тригонометричні, показникові, логарифмічні і т.д. функції, то рівняння (1) називають трансцендентним (нелінійним).
Знайти точні значення коренів заданого рівняння можна лише для найпростіших функцій f(x); алгебраїчних многочленiв не вище четвертого степеня, деяких многочленiв степеня n5 і деяких трансцендентних функцій.
Універсальних методів для знаходження точних значень коренів алгебраїчних рівнянь степеня n5 і трансцендентних рівнянь не існує. Тому важливого значення набувають наближені методи знаходження коренів з достатньою для практики точністю. Задача знаходження коренів рівняння (1) вважається розв'язаною, якщо корені обчислені із наперед заданою точністю.
Знаходження наближених коренів рівняння (1) складається з двох етапів:
1) відокремлення коренів, тобто знаходження досить малих відрізків, на кожному з яких міститься один і тільки один корінь рівняння;
2) обчислення коренів з наперед заданою точністю.
До всіх представлених методів розроблені програмні реалізації на мові С. Для порівняння, завдання виконані також за допомогою пакету MathCad.
Методи хорд, дотичних, хорд-дотичних
В даному розділі реалізовані такі методи
метод хорд
метод дотичних
метод хорд-дотичних
Ці методи полягають у наближеній заміні на досить малому відрізку функції f(x) лінійною функцією, яка відповідає хорді, проведеній через дві задані точки (метод хорд), дотичною, проведеною в заданій точці (метод дотичних), або ж ці методи комбінуються (метод хорд-дотичних).
Метод хорд
Якщо на інтервалі [a;b] неперервна функція F(x) задовольняє умову F(a)*F(b)<0, то корінь рівняння F(x)=0 наближено знаходиться за рекурентною формулою
xn=xn-1- (2)
Похибка обчислюється за формулою
|xn-xn-1|< (3)
Рис.1
Метод дотичних
Корінь рівняння F(x)=0 обчислюється за ітераційною формулою
xk+1=xk- (4)
Графічна ілюстрація приведена на рис.2.
рис.2
Метод хорд-дотичних
При використанні комбінованого методу хорд-дотичних уточнення інтервалу проводиться як методом хорд, так і методом дотичних. В залежності від знаку функції на кінцях уточненого інтервалу відбувається вибір наступного інтервалу. Графічна ілюстрація методу приведена на рис.3.
рис.3
Завдання
Використовуючи методи хорд, дотичних, хорд-дотичних знайти корінь функції f(x)=x3-2x2+x-3 на інтервалі [2,1;2,2].
Програмна реалізація методів
Метод хорд
Текст програми
#include<stdio.h>
#include<math.h>
/*метод хорд*/
float fun(float x) {
return(x*x*x-2*x*x+x-3);
};
main() {
float e=0.00001,a=2.1,b=2.2,a1,r;
FILE *res;
res=fopen("resch.txt","w");
fprintf(res,"eps=%5.6f\n",e);
do {
a1=a-fun(a)*(b-a)/(fun(b)-fun(a));
if (fun(a1)<0) a=a1;
else b=a1;
r=fabs(fun(a));
fprintf(res,"a=%5.3f b=%5.3f f(a)=%5.3f\n",a,b,r);
} while(r>=e);
fprintf(res,"корiнь=%5.3f f(x)=%5.3f\n",a,r);
fclose(res);
}
Результат
eps=0.00001
a=2.173 b=2.200 f(a)=0.009
a=2.175 b=2.200 f(a)=0.000
a=2.175 b=2.200 f(a)=0.000
корiнь=2.175 f(x)=0.000
Метод дотичних
Текст програми
#include<stdio.h>
#include<math.h>
/*метод дотичних*/
float fun(float x) {
return(x*x*x-2*x*x+x-3);
};
float funp(float x) {
return(3*x*x-4*x+1);
};
main() {
float e=0.00001,a=2.1,b=2.2,b1,r;
FILE *res;
res=fopen("resdot.txt","w");
fprintf(res,"eps=%5.6f\n",e);
do {
b1=b-fun(b)/funp(b);
if (fun(b1)>0) b=b1;
else a=b1;
r=fabs(fun(b));
fprintf(res,"a=%5.3f b=%5.3f f(b)=%5.3f\n",a,b,r);
} while(r>=e);
fprintf(res,"корiнь=%5.3f f(x)=%5.3f\n",b,fun(b));
fclose(res);
}
Результат
eps=0.00001
a=2.100 b=2.175 f(b)=0.003
a=2.100 b=2.175 f(b)=0.000
корiнь=2.175 f(x)=0.000
Метод хорд-дотичних
Текст програми
#include<stdio.h>
#include<math.h>
/*метод хорд-дотичних*/
float fun(float x) {
return(x*x*x-2*x*x+x-3);
};
float funp(float x) {
return(3*x*x-4*x+1);
};
main() {
float e=0.00001,a=2.1,b=2.2,c;
FILE *res;
res=fopen("reschdot.txt","w");
do {
a=a-fun(a)*(b-a)/(fun(b)-fun(a));
b-=fun(b)/funp(b);
fprintf(res,"a=%5.3f b=%5.3f b-a=%5.3f\n",a,b,b-a);
} while(b-a>=e);
c=(b+a)/2;
fprintf(res,"корiнь=%5.3f f(x)=%5.3f\n",c,fun(c));
fclose(res);
}
Результат
a=2.173 b=2.175 b-a=0.002
a=2.175 b=2.175 b-a=0.000
корiнь=2.175 f(x)=-0.000
Метод ітерацій
Для застосування методу ітерацій вихідне рівняння слід записати у формі
x=(x) (4).
Нехай виділений інтервал <a;b> ізоляції кореня цього рівняння, і x0 - довільна точка цього інтервалу (нульове наближення). Для одержання наступного наближення в праву частину рівняння (4) замість x підставляємо x0 і т.д. Ітераційний процес збігається, якщо відображення (x) - стягуюче.
Завдання
За допомогою методу ітерацій знайти корінь функції f(x)= при нульовому наближенні 2,1.
Програмна реалізація методу
Текст програми
#include<stdio.h>
#include<math.h>
/*метод iтерацiй*/
float fun(float x) {
return(1/(1+exp(x)));
};
main() {
float e=0.00001,a=2.1,r,x;
FILE *res;
x=a;
res=fopen("resit.txt","w");
fprintf(res,"eps=%5.6f\n",e);
do {
r=x;
x=fun(x);
r=fabs(r-x);
} while(r>=e);
fprintf(res,"корiнь=%5.3f \n",x);
fclose(res);
}
Результат
eps=0.00001
корiнь=0.401
Метод половинного поділу
Метод поділу відрізка пополам (метод дихотомії) застосовний для уточнення кореня рівняння f(x) з наперед заданою точністю, якщо функція f(x) задовольняє такі умови: неперервність і диференційованість на відрізку <a;b>, набуває на кінцях цього відрізка значень різних знаків, а похідна f'(x) зберігає сталий знак всередині відрізка <a;b>. Уточнення інтервалу відбувається методом половинного поділу.
Завдання
За допомогою методу половинного поділу знайти корінь функції f(x)=x3-x2+1 на відрізку [2;3].
Програмна реалізація методу
Текст програми
/*пошук коpенiв методом половинного подiлу*/
#include <stdio.h>
float f1(float x) {
return(x*x*x-x*x+1);
}
int sign(float x) {
return(x>0?1:x?-1:0);
}
main() {
float a,b,p,p1,eps,x,d;
FILE *stream;
stream=fopen("der\\kor.txt","w");
printf("Введiть межi a,b\n");
scanf("%f",&a);
scanf("%f",&b);
fprintf(stream,"Для функцii f=x^3-x^2+1\n");
fprintf(stream,"В межах %5.2f %5.2f\n",a,b);
printf("Введiть точнiсть обчислень e=");
scanf("%f",&eps);
p=(b-a)/2;
p1=(a+b)/2;
for (x=p1,d=p;d>eps;d/=2) x-=sign(f1(x))*d;
fprintf(stream,"f(%5.3f)=0 з точнiстю %5.3f\n",x,eps);
fclose(stream);
}
Результат
Для функцii f=x^3-x^2+1
В межах -2.00 3.00
f(-0.756)=0 з точнiстю 0.01
МЕТОДИ ЛІНІЙНОЇ АЛГЕБРИ
Знаходження визначника матриці методом алгебраїчних доповнень
В даному розділі представлено метод знаходження визначника матриці за допомогою мінорів.
Введемо такі позначення: якщо aij - елемент визначника d, то через Mij позначимо додатковий мінор, чи, коротше, мінор цього елемента, тобто мінор (n-1)-го порядку, який одержуємо після викреслення з визначника i-го рядка та j-го стовпця. Далі, через Aij позначимо алгебраїчне доповнення елемента aij, тобто
Aij=(-1)i+jMij
Визначник d рівний сумі добутків всіх елементів довільного його рядка на їх алгебраїчні доповнення.
d=ai1Ai1+ai2Ai2+...+ainAin (1)
Заміняючи у розкладі (1) алгебраїчні доповнення відповідними мінорами зі знаками плюс чи мінус, зведемо обчислення визначника n-го порядку до обчислення кількох визначників (n-1)-го порядку. При програмуванні можна використати рекурсію.
Представлена програма, написана на мові С, здійснює обчислення визначника. Представлено результат роботи програми. Саме обчислення визначника оформлено у вигляді функції, яка здійснює рекурсію - викликає сама себе.
Текст програми
#include <stdio.h>
#include <math.h>
/*пiдпрограма обчислення визначника*/
/*методом алгебраiчних доповнень*/
float det(float a[10][10],int n) {
int j,k,l,k1,l1,zn,n1;
float b[10][10],deter,d;
if (n==1) deter=a[0][0];
if (n==2) deter=a[0][0]*a[1][1]-a[0][1]*a[1][0];
if (n>2) {
deter=0;
for (j=0;j<n;j++) {
if (fmod(j,2)==0) zn=1;
else zn=-1;
k=1;l=0;n1=n-1;
for (k1=0;k1<n1;k1++,k++) {
l=0;
for (l1=0;l1<n1;l1++,l++) {
if (l1==j) l++;
b[k1][l1]=a[k][l];
};
};
d=det(b,n1);
deter+=a[0][j]*zn*d;
};
};
return(deter);
}
main () {
float a,mas[10][10];
int i,j,p;
FILE *res;
res=fopen("deter.txt","w");
printf("Введiть порядок визначника ");
scanf("%d",&p);
/*ввiд коефiцiентiв рiвняннь*/
for(i=0;i<p;i++) {
for(j=0;j<p;j++) {
scanf("%f",&a);
mas[i][j]=a;
};
};
fprintf(res,"Матриця\n");
for(i=0;i<p;i++) {
for(j=0;j<p;j++) fprintf(res," %5.2f",mas[i][j]);
fprintf(res,"\n");
};
a=det(mas,p);
fprintf(res,"визначник матрицi det=%5.3f\n ",a);
fclose(res);
}
Результат роботи програми
Матриця
3.00 6.00 4.00 1.00
7.00 0.00 4.00 3.00
2.00 5.00 7.00 9.00
3.00 7.00 4.00 2.00
визначник матрицi det=-308.000
Розв'язок системи рівнянь методом Крамера
Один з методів розв'язку системи лінійних рівнянь - правило Крамера. Використовуючи позначення попереднього розділу, вони формулюються так: якщо система лінійних рівнянь сумісна, то вона має єдиний розв'язок
a1=
де dj - визначник, одержаний із визначника системи d заміною j-го стовпця стовпцем із вільних членів системи рівнянь.
Представлена програма на мові С реалізує метод Крамера. Визначники обчислюються за тим же правилом, що й у попередньому розділі - з використанням рекурсії. Для порівняння, та ж система рівнянь розв'язана методом Крамера за допомогою пакету MathCad. Результати співпадають.
Текст програми
/* Програма kra.c */
/* Програма призначена для розв'язку системи рiвнянь */
/* методом Крамера. */
#include <stdio.h>
#include <math.h>
/*пiдпрограма обчислення визначника*/
float det(float a[10][10],int n) {
int j,k,l,k1,l1,zn,n1;
float b[10][10],deter,d;
if (n==1) deter=a[0][0];
if (n==2) {
deter=a[0][0]*a[1][1]-a[0][1]*a[1][0];
};
if (n>2) {
deter=0;
for(j=0;j<n;j++) {
if (fmod(j,2)==0) zn=1;
else zn=-1;
k=1;l=0;n1=n-1;
for (k1=0;k1<n1;k1++,k++) {
l=0;
for (l1=0;l1<n1;l1++,l++) {
if (l1==j) l++;
b[k1][l1]=a[k][l];
/*printf("b%d%d=%5.2f a%d%d=%5.2f\n",k1,l1,b[k1][l1],k,l,a[k][l]);*/
};
};
d=det(b,n1);
deter+=a[0][j]*zn*d;
/*printf("n1=%d deter=%.2f det=%.2f\n",n1,deter,d);*/
};
};
return(deter);
}
main () {
float a,c,mas[10][10],mas1[10][10],b[10],x[10],r[10],s;
int i,j,p,k,l;
char d[10]="введiть a";
char e[10]="введiть b";
FILE *potik;
printf("Введiть порядок системи ");
scanf("%d",&p);
/*ввiд коефiцiентiв рiвняннь*/
for (i=0;i<p;i++) {
for (j=0;j<p;j++) {
printf("%s %d %d\n",d,i+1,j+1);
scanf("%f",&a);
mas[i][j]=a;
};
};
/*ввiд вiльних членiв*/
for (i=0;i<p;i++) {
printf("%s %d\n",e,i+1);
scanf("%f",&a);
b[i]=a;
};
potik=fopen("resk.txt","w");
fprintf(potik," Розв'язок системи\n");
fprintf(potik," за формулами Крамера\n");
fprintf(potik," матриця системи\n");
for (i=0;i<p;i++) {
for (j=0;j<p;j++) fprintf(potik,"%.2f ",mas[i][j]);
fprintf(potik,"=%.2f\n",b[i]);
};
a=det(mas,p);
fprintf(potik,"визначник системи det=%5.3f\n ",a);
/*перевiрка невиродженостi системи*/
if (a==0)
fputs("система невизначена",potik);
else {
fprintf(potik," результати\n");
/*вивiд результатiв*/
for (i=0;i<p;i++) {
for (k=0;k<p;k++) for(l=0;l<p;l++) mas1[k][l]=mas[k][l];
for(j=0;j<p;j++) mas1[j][i]=b[j];
c=det(mas1,p);
x[i]=c/a;
fprintf(potik,"x%d= %.3f\n",i,x[i]);
};
fprintf(potik," невязки\n");
for (i=0;i<p;i++) {
s=0;
for (j=0;j<p;j++){
s+=mas[i][j]*x[j];
};
r[i]=s-b[i];
fprintf(potik,"r[%d]=%5.3f\n",i,r[i]);
};
};
fclose(potik);
}
Розв'язок системи
за формулами Крамера
матриця системи
2.00 1.00 1.00 =2.00
4.00 3.00 2.00 =2.00
1.00 4.00 5.00 =1.00
результати
x0= 1.222
x1= -2.000
x2= 1.556
невязки
r[0]=0.000
r[1]=0.000
r[2]=0.000
Розв'язок системи лінійних рівнянь методом Крамера
за допомогою пакету MathCad
Розв'язок системи рівнянь методом Гауса
Розв'язок системи лінійних рівнянь за правилами Крамера громіздкий і довготривалий, навіть з використанням ПЕОМ. Найпростішим методом розв'язування систем лінійних алгебраїчних рівнянь є метод послідовного виключення змінних або метод Гауса. На першому етапі вихідну систему рівнянь зводять до рівносильної їй системи трикутної форми. Цей процес перетворення називають прямим ходом. На другому етапі, який називають зворотним ходом, знаходять розв'язок лінійної системи рівнянь трикутної форми.
Представлена програма на мові Pascal (версія TurboPascal6 і вище) здійснює реалізацію одної з модифікацій методу Гауса з вибором головного елемента. Сам метод Гауса оформлений як окрема процедура. Представлено результат роботи програми з тією ж системою, яка була розв'язана за правилом Крамера, а також розв'язок тієї ж системи засобами пакету MathCad. Результати співпадають.
Текст програми
program gauss(input,output);
{розв'язування системи рiвнянь методом гауса'}
TYPE
dvmr = array[0..10,0..10] of real;
ovmr = array[0..10] of real;
var i,j,n,m:integer;
n1,m1,l1:integer;
aa,a,p,g:dvmr;
x,y,riz:ovmr;
b,ii:integer;
ab:real;
fr:text;
Procedure GELG(A:dvmr; Var B:ovmr; n:integer);
Var i,j,k,l : integer;
r : real;
Begin
{ ************************************************}
{ Пiдпрограма розв`язування системи лiнiйних }
{ алгебраїчних рiвнянь методом Гауса }
{ з вибором головного елемента }
l:=0;
for i:=1 to N do begin
k:=i; r:=Abs(A[i,i]);
for j:=i+1 to N do
if Abs(A[j,i])>r then begin
k:=j; r:=Abs(A[j,i]);
end;
if r<>0 then begin
if k<>i then begin
r:=B[k]; B[k]:=B[i]; B[i]:=r;
for j:=i to N do begin
r:=A[k,j];
A[k,j]:=A[i,j];
A[i,j]:=r;
end;
end;
r:=A[i,i]; B[i]:=B[i]/r;
for j:=i to N do A[i,j]:=A[i,j]/r;
for k:=i+1 to N do begin
r:=A[k,i]; B[k]:=B[k]-r*B[i];
for j:=i to N do A[k,j]:=A[k,j]-r*A[i,j];
end;
end
else begin l:=1; i:=N+1; end;
end;
if l<>1 then begin
for i:=N-1 downto 1 do begin
for j:=i+1 to N do B[i]:=B[i]-A[i,j]*B[j];
end;
end;
End;
{процедура вводу матрицi }
procedure vvidmat(var a:dvmr;n,m:integer);
var i,j:integer;
begin
for i:=1 to n do
for j:=1 to m do
readln(a[i,j]);
end;
{процедура вводу вектора }
procedure vvidvek(var x:ovmr;n:integer);
var i:integer;
begin
for i:=1 to n do
readln(x[i]);
end;
{процедура виводу матрицi на друк}
procedure drukmat(a:dvmr;n,m:integer);
var i,j:integer;
begin
for i:=1 to n do
begin
for j:=1 to m do
write(fr,a[i,j]:5:2,' ');
writeln(fr);
end;
end;
procedure drukvek(x:ovmr;n:integer);
var i:integer;
begin
for i:=1 to n do
writeln(fr,x[i]:5:3);
end;
begin
assign(fr,'gres.txt');
rewrite(fr);
{побудова матрицi системи}
writeln('введiть порядок системи');
readln(m);
writeln('ввiд с-ми');
vvidmat(aa,m,m);
writeln('ввiд вiльних членiв');
vvidvek(x,m);
rivvek(x,y,m);
writeln(fr,'матриця системи');
drukmat(aa,m,m);
writeln(fr,'вiльнi члени');
drukvek(x,m);
gelg(aa,x,m);
{знаходження i друк розв"язкiв}
writeln(fr,'розв"язок');
drukvek(x,m);
{знаходження нев'язок}
writeln(fr,'нев"язки');
for i:=1 to m do
begin
riz[i]:=0;
for j:=1 to m do
riz[i]:=riz[i]+aa[i,j]*x[j];
riz[i]:=riz[i]-y[i];
writeln(fr,'riz[',i:1,']=',riz[i]:2:2);
end;
close(fr);
end.
Результат роботи програми
матриця системи
2.00 1.00 1.00
4.00 3.00 2.00
1.00 4.00 5.00
вiльнi члени
2.00
2.00
1.00
розв"язок
1.222
-2.000
1.556
нев"язки
riz[1]=0.00
riz[2]=0.00
riz[3]=0.00
Розв'язок системи рівнянь за допомогою пакету MathCad
Знаходження визначника матриці методом Гауса
В результаті перетворень Гауса система (або матриця) зводиться до форми, яка дає змогу легко обчислити визначник матриці.
Представлена програма на мові Pascal (версія TurboPascal6 і вище) здійснює перетворення матриці за схемою Гауса та знаходження визначника. Для контролю визначник тієї ж матриці обчислено за допомогою пакету MathCad. Результати співпадають.
Текст програми
program gauss6(input,output);
{знаходження визначника методом Гауса}
TYPE
dvmr = array[0..10,0..10] of real;
ovmr = array[0..10] of real;
var i,j,n,m:integer;
n1,m1,l1:integer;
aa,a,p,g:dvmr;
x,y,riz:ovmr;
b,ii:integer;
d,ab:real;
fr:text;
Function Det(A:dvmr; N:integer):real;
Var i,j,k : integer;
d,y,w : real;
Begin
{**************************************************}
{* Знаходження визначника матрицi методом Гауса *}
d:=1;
for i:=1 to N do begin
k:=i; y:=A[i,i];
for j:=i+1 to N do begin
w:=A[j,i];
if Abs(w)>Abs(y) then begin
k:=j; y:=w;
end;
end;
d:=d*y;
if d=0 then begin det:=0; Exit; end;
if i<>k then begin
d:=-d;
for j:=i to N do begin
w:=A[k,j];
A[k,j]:=A[i,j];
A[i,j]:=w;
end;
end;
for j:=i+1 to N do begin
w:=A[j,i]/y;
for k:=i+1 to N do A[j,k]:=A[j,k]-w*A[i,k];
end;
end;
Det:=d;
End;
{процедура вводу матрицi }
procedure vvidmat(var a:dvmr;n,m:integer);
var i,j:integer;
begin
for i:=1 to n do
for j:=1 to m do
readln(a[i,j]);
end;
{процедура виводу матрицi на друк}
procedure drukmat(a:dvmr;n,m:integer);
var i,j:integer;
begin
for i:=1 to n do
begin
for j:=1 to m do
write(fr,a[i,j]:5:2,' ');
writeln(fr);
end;
end;
begin
assign(fr,'det.txt');
rewrite(fr);
{побудова матрицi }
writeln('введiть порядок матрицi');
readln(m);
writeln('ввiд с-ми');
vvidmat(aa,m,m);
writeln(fr,'матриця ');
drukmat(aa,m,m);
d:=det(aa,m);
{друк визначника}
writeln(fr,'визначник=',d:5:2);
writeln(fr,'======================================================');
close(fr);
end.
Результат роботи програми
матриця
3.00 6.00 4.00 1.00
7.00 0.00 4.00 3.00
2.00 5.00 7.00 9.00
3.00 7.00 4.00 2.00
визначник=-308.00
Знаходження визначника матриці за допомогою пакету MathCad
Знаходження оберненої матриці за допомогою методу Гауса
Перетворивши матрицю за схемою Гауса, легко знайти обернену матрицю.
Представлена програма на мові Pascal здійснює перетворення Гауса (зведення матриці до трикутної форми), обчислює визначник матриці та знаходить обернену матрицю. Представлено результат роботи програми та розв'язок цієї ж задачі засобами пакету MathCad. Результати співпадають.
Текст програми
program gauss(input,output);
{знаходження визначника та оберненоi матрицi методом Гауса}
TYPE
dvmr = array[0..10,0..10] of real;
ovmr = array[0..10] of real;
var i,j,n,m:integer;
n1,m1,l1:integer;
aa,a,p,g:dvmr;
x,y,riz:ovmr;
b,ii:integer;
d,ab:real;
fr:text;
Procedure Inv_Gauss(Var A:dvmr; N:integer;Var D:real);
Var i,j,k : integer;
y,w,e : real;
z : array[1..10] of integer;
Begin
{**************************************************}
{ Знаходження оберненоi матрицi методом Гауса *}
D:=1; e:=0.00001;
for j:=1 to N do z[j]:=j;
for i:=1 to N do begin
k:=i; y:=A[i,i];
for j:=i+1 to N do begin
w:=A[i,j];
if Abs(w)>Abs(y) then begin k:=j; y:=w; end;
end;
D:=y*D;
if Abs(y)<e then begin
writeln('Матриця вироджена');
Exit;
end;
for j:=1 to N do begin
A[0,j]:=A[j,k]; A[j,k]:=A[j,i];
A[j,i]:=-A[0,j]/y;
A[i,j]:=A[i,j]/y;
A[j,0]:=A[i,j];
end;
A[i,i]:=1/y; j:=z[i]; z[i]:=z[k]; z[k]:=j;
for k:=1 to N do
if k<>i then
for j:=1 to N do
if j<>i then A[k,j]:=A[k,j]-A[j,0]*A[0,k];
end;
for i:=1 to N do begin
k:=z[i];
if k<>i then begin
for j:=1 to N do begin
w:=A[i,j];
A[i,j]:=A[k,j];
A[k,j]:=w;
end;
j:=z[i]; z[i]:=z[k]; z[k]:=j; D:=-D;
end;
end;
End;
{процедура вводу матрицi }
procedure vvidmat(var a:dvmr;n,m:integer);
var i,j:integer;
begin
for i:=1 to n do
for j:=1 to m do
readln(a[i,j]);
end;
{процедура виводу матрицi на друк}
procedure drukmat(a:dvmr;n,m:integer);
var i,j:integer;
begin
for i:=1 to n do
begin
for j:=1 to m do
write(fr,a[i,j]:5:2,' ');
writeln(fr);
end;
end;
begin
assign(fr,'obmat.txt');
rewrite(fr);
{побудова матрицi }
writeln('введiть порядок матрицi');
readln(m);
writeln('ввiд матрицi');
vvidmat(aa,m,m);
writeln(fr,'матриця ');
drukmat(aa,m,m);
inv_gauss(aa,m,d);
{друк визначника}
writeln(fr,'визначник=',d:5:2);
writeln(fr,'======================================================');
writeln(fr,'обернена матриця');
drukmat(aa,m,m);
close(fr);
end.
Результат роботи програми
матриця
3.00 1.00 5.00 2.00
6.00 -1.00 3.00 2.00
1.00 7.00 4.00 3.00
2.00 0.00 7.00 8.00
визначник=747.00
======================================================
обернена матриця
-0.108 0.229 0.048 -0.048
-0.078 0.028 0.158 -0.047
0.415 -0.185 -0.086 -0.025
-0.336 0.104 0.063 0.159
Знаходження оберненої матриці за допомогою пакету MathCad
Розв'язок системи лінійних рівнянь за допомогою методу Зейделя
Застосування методу Гауса для розв'язування системи лінійних рівнянь з великою кількістю невідомих досить громіздке. Крім того, кількість невідомих може бути така велика, що коефіцієнти системи не завжди можна розмістити у оперативній пам'яті ПЕОМ. Тоді застосовувати для її розв'язування метод Гауса взагалі не можна. У цих випадках розв'язують систему ітераційними методами. Метод Зейделя належить саме до ітераційних методів.
Представлена програма на мові С здійснює реалізацію методу Зейделя.
Текст програми
#include <stdio.h>
#include <math.h>
main()
{
FILE *stream;
float e=0.001,a[11][11],x[11],y[11],s,r[11],t[11][11];
int i,j,n,m;
stream=fopen("zej.txt","w");
puts("введiть к-сть рiвнянь");
scanf("%d",&m);
m--;
for (i=0;i<=m;i++)
for (j=0;j<=m+1;j++)
{printf("a %d %d=",i,j);
scanf("%f",&a[i][j]);
t[i][j]=a[i][j];};
fprintf(stream," Розв'язок системи\n");
fprintf(stream," методом Зейделя\n");
fprintf(stream," матриця системи\n");
for (i=0;i<=m;i++)
{
for (j=0;j<=m+1;j++)
fprintf(stream," %5.2f",a[i][j]);
fprintf(stream,"\n");
};
/*присвоення початкових значень*/
for(i=0;i<=m;i++)
x[i]=a[i][m+1];
/*цикл iтерацiй*/
do {
s=0;
for(n=0;n<=m;n++)
{
y[n]=a[n][m+1];
for(j=0;j<=m;j++)
{if (j==n) continue;
y[n]-=a[n][j]*x[j];};
y[n]/=a[n][n];
s+=fabs(x[n]-y[n]);
x[n]=y[n];
};}
while (s>=e);
fprintf(stream," Розв'язок\n");
for (i=0;i<=m;i++)
fprintf(stream,"x %d=%5.4f\n",i,x[i]);
fprintf(stream,"невязки\n");
for(i=0;i<=m;i++)
{s=0;
for(j=0;j<=m;j++)
{
s+=t[i][j]*x[j]; };
r[i]=s-t[i][m+1];
fprintf(stream,"r[%d]=%5.3f\n",i,r[i]);
};
fclose(stream);
}
Результат роботи програми
Розв'язок системи
методом Зейделя
матриця системи
4.00 0.10 0.20 5.00
-0.20 7.00 0.60 5.00
-0.20 -0.30 9.00 12.00
Розв'язок
x 0=1.1653
x 1=0.6293
x 2=1.3802
невязки
r[0]=0.000
r[1]=0.000
r[2]=0.000
ІНТЕРПОЛЮВАННЯ ФУНКЦІЙ
Нехай на відрізку [a;b] визначено певний клас функцій {P(x)}, наприклад, клас алгебраїчних многочленiв, а в точках x0,x1,...,xn цього проміжку задано значення деякої функції y=f(x): y0=f(x0), y1=f(x1), ..., yn=f(xn). Наближену заміну функції f на відрізку [a;b] однією з функцій P(x) цього класу так, щоб функція P(x) в точках x0,x1,...,xn набувала тих самих значень, що й функція f, називають інтерполюванням або інтерполяцією. Точки x0, x1, ... ,xn називають вузлами інтерполювання, функцію P(x) - інтерполюючою функцією, а формулу f(x)P(x), за допомогою якої обчислюють значення функції f у проміжку [a;b], - інтерполяційною формулою.
Якщо функція P(x) належить до класу алгебраїчних многочленiв, то інтерполювання називається параболічним. Параболічне інтерполювання найзручніше, оскільки многочлени, які прості за формою і не мають особливих точок, можуть набувати довільних значень, їх легко обчислювати, диференціювати та інтегрувати.
Сформулюємо задачу параболічного інтерполювання: в n+1 різних точках x0, x1, ... ,xn задано значення функції f: y0=f(x0), y1=f(x1), ..., yn=f(xn) і треба побудувати многочлен
Pn(x)=a0xn+a1xn-1+...+an-1x+an
степеня n, який задовольняв би умови:
Pn(xi)=yi (i=0,1, . . . , n).
Задача має єдиний розв'язок. Многочлен Pn(x) називають інтерполяційним многочленом. Інтерполяційний многочлен єдиний, проте можливі різні форми його запису.
Інтерполяційний многочлен будують тоді, коли:
функцію задано таблично для деяких значень аргументу, а треба знайти її значення для значень аргументу, яких в таблиці нема.
функцію задано графічно, а треба знайти її наближений аналітичний вираз.
функцію задано аналітично, але її вираз досить складний і незручний для виконання різних математичних операцій.
При написанні даної роботи розглядалася перша задача - чисельної інтерполяції.
Інтерполяційний многочлен Лагранжа
Інтерполяційний многочлен Лагранжа має такий вираз:
Ln(x)=
Многочлен Лагранжа зручно будувати у випадку рівновіддалених вузлів.
Завдання. Задане табличне представлення функції:
|
x |
3 |
7 |
11 |
15 |
19 |
|
|
y |
4 |
10 |
22 |
26 |
23 |
Методом інтерполяції Лагранжа знайти значення функції при x=13.
Програмна реалізація здійснена на мові С.
Текст програми
#include<stdio.h>
#include<math.h>
/* iнтерполяцiя ф-ii за формулою Лагранжа при N вузлах*/
main()
{
float x[10],y[10],l,p,z,t;
int i,n,j;
FILE *stream;
stream=fopen("lagr.txt","w");
fprintf(stream,"Iнтерполяцiя за Лагранжем\n");
printf("Введiть к-сть вузлiв ");
scanf("%d",&n);
printf("Введiть значення x ");
for(i=0;i<n;i++) scanf("%f",&x[i]);
printf("Введiть значення y ");
for(i=0;i<n;i++) scanf("%f",&y[i]);
printf("Введiть значення т-ки iнтерполяцii");
scanf("%f",&t);
for (i=0;i<n;i++) fprintf(stream,"%5.3f %5.3f\n",x[i],y[i]);
l=0;
for(i=0;i<n;i++)
{
p=1;
for(j=0;j<n;j++)
{
if(i==j) z=t-x[j];
else
z=x[i]-x[j];
p*=(t-x[j])/z;
};
l+=p*y[i];
};
fprintf(stream,"\nзначення полiнома Лагранжа при x=%5.3f = %5.4f\n",t,l);
fclose(stream);
}
Результат роботи програми
Iнтерполяцiя за Лагранжем
x y
3.000 4.000
7.000 10.000
11.000 22.000
15.000 26.000
19.000 33.000
значення полiнома Лагранжа при x=13.000 = 24.8984
Побудова інтерполяційного многочлена Лагранжа
Дещо більшу складність порівняно з числовою інтерполяцією має задача побудови самого многочлена. Ця задача розв'язувалася засобами Pascal (версія TurboPascal 6.0 і вище). При написанні далі приведеної програми використовувався модуль операцій з многочленом bibl.tpu. Процедури і функції, які входять до цього модуля, дають змогу здійснювати операції додавання, віднімання, множення, ділення з остачею і ще деякі в кільці многочленів.
program lagr;
{$M 65520,0,655360}
{побудова многочлена Лагранжа}
Uses Crt,bibl;
{початок програми}
var i,j,k,n,m:integer;
s,p,q,p1:poli;
t,u,w:real;
x,y:array[1..20] of real;
begin
{створення кiльцевого нуля zero i кiльцевоi одиницi od}
zerod;
assign(fi,'lagr.txt');
rewrite(fi);
{ввiд вузлiв}
writeln('Введiть число вузлiв ');
readln(n);
for i:=1 to n do begin
writeln('Введiть x[',i,'] y[',i,']');
readln(x[i],y[i]); end;
writeln('Введiть точку iнтерполяцii '); readln(t);
writeln(' x y');
for i:=1 to n do writeln(x[i]:5:2,' ',y[i]:5:2);
writeln(fi,' x y');
for i:=1 to n do writeln(fi,x[i]:5:2,' ',y[i]:5:2);
writeln('Точка iнтерполяцii ',t:5:3);
writeln(fi,'Точка iнтерполяцii ',t:5:3);
s:=zero;
for i:=1 to n do
begin
p:=od;u:=1;
for k:=1 to n do
begin
if (k<>i) then
begin
q:=zero;q[1]:=1;q[0]:=-x[k];
dobutok(p,q,p); u:=u*(x[i]-x[k]);
end; end;
dobchy(p,y[i]/u,p);
suma(s,p,s);end;
writeln('Многочлен Лагранжа '); writeln(fi,'Многочлен Лагранжа ');
vyvid(s); fvyvid(s);
writeln; writeln(fi);
w:=znach(s,t);
writeln('Значення в точцi iнтерполяцii=',w:5:3);
writeln(fi,'Значення в точцi iнтерполяцii=',w:5:3);
close(fi);
end.
Результат роботи програми
x y
3.00 4.00
7.00 10.00
11.00 22.00
15.00 26.00
19.00 33.00
Точка iнтерполяцii 13.000
Многочлен Лагранжа
0.004x^ 4-0.183x^ 3+2.768x^ 2-14.087x+25.958
Значення в точцi iнтерполяцii=24.898
Інтерполяційний многочлен Ньютона
Інтерполяційний многочлен Ньютона будується з використанням скінченних різниць . Його перевага порівняно з інтерполяційним многочленом Лагранжа полягає у простішій модифікації при добавленні ще одного вузла.
Вигляд інтерполяційного многочлену Ньютона:
Завдання. Задане табличне представлення функції:
|
X |
3 |
7 |
11 |
15 |
19 |
|
|
Y |
4 |
10 |
22 |
26 |
23 |
Методом інтерполяції Нютона знайти значення функції при x=13.
Програмна реалізація
Здійснена на мові С.
Текст програми
/*iнтерполяцiя за методом Ньютона*/
#include<stdio.h>
main()
{
FILE *stream;
int i,j,n;
float x[10],y[10],d[10][10],xx,h,nh,p,s,fact;
printf("Введiть к-сть вузлiв iнтерполяцii");
scanf("%d",&n);
stream=fopen("ne.txt","w");
fprintf(stream,"Iнтерполяцiя за Ньютоном\n");
fprintf(stream,"x ");
for (i=0;i<n;i++)
{puts("ввiд x");
scanf("%f",&x[i]);
fprintf(stream,"%5.2f ",x[i]);
};
fprintf(stream,"\ny ");
for (i=0;i<n;i++)
{puts("ввiд y");
scanf("%f",&y[i]);
fprintf(stream,"%1.3f ",y[i]);
};
printf("Введiть значення точки iнтерполяцii ");
scanf("%f",&xx);
for(i=0;i<=n-1;i++)
d[i][0]=y[i+1]-y[i];
for(j=1;j<n;j++)
for(i=0;i<n;i++)
d[i][j]=d[i+1][j-1]-d[i][j-1];
fprintf(stream,"\nx=%5.3f",xx);
s=y[0];
p=1.0;
h=x[1]-x[0];
nh=1.0;
fact=1.0;
for(i=0;i<n-1;i++)
{ fact*=(i+1);
nh*=h;
p*=xx-x[i];
s+=d[0][i]*p/(fact*nh);
};
fprintf(stream,"\ny=%5.4f",s);
printf("y=%5.4f",s);
fclose(stream);
}
Результат роботи програми
Iнтерполяцiя за Ньютоном
x 3.00 7.00 11.00 15.00 19.00
y 4.000 10.000 22.000 26.000 33.000
x=13.000
y=24.8984
Результат спіпав з результатом, одержаним за допомогою многочлена Лагранжа.
Побудова інтерполяційного многочлена Ньютона
Побудова інтерполяційного многочлена Ньютона здійснювалася засобами Pascal з використанням модуля роботи в кільці многочленів bibl.tpu, вихідний текст якого приведений у розділі “Побудова многочлена Лагранжа”. Далі приведений текст програми побудови многочлена Ньютона.
program newt;
{$M 65520,0,655360}
{побудова многочлена Ньютона}
Uses Crt,bibl;
{початок програми}
var i,j,k,n,m:integer;
s,p,q,p1:poli;
fact,h,nh,t,w:real;
x,y:array[1..20] of real;
d:array[1..20,1..20] of real;
begin
{створення кiльцевого нуля zero i кiльцевоi одиницi od}
zerod;
assign(fi,'newt.txt');
rewrite(fi);
{ввiд вузлiв}
writeln('Введiть число вузлiв ');
readln(n);
for i:=1 to n do begin
writeln('Введiть x[',i,'] y[',i,']');
readln(x[i],y[i]);
end;
writeln('Введiть точку iнтерполяцii ');
readln(t);
writeln(' x y');
for i:=1 to n do writeln(x[i]:5:2,' ',y[i]:5:2);
writeln(fi,' x y');
for i:=1 to n do writeln(fi,x[i]:5:2,' ',y[i]:5:2);
writeln('Точка iнтерполяцii ',t:5:3);
writeln(fi,'Точка iнтерполяцii ',t:5:3);
s:=zero;
s[0]:=y[1];
h:=x[2]-x[1];
for i:=1 to n do d[i,1]:=y[i+1]-y[i];
for j:=2 to n do
for i:=1 to n do
d[i,j]:=d[i+1,j-1]-d[i,j-1];
fact:=1;nh:=1;
p:=od;
for i:=1 to n-1 do
begin
fact:=fact*i;
nh:=nh*h;
q:=zero;
q[1]:=1;q[0]:=-x[i];
dobutok(p,q,p);
dobchy(p,d[1,i]/fact/nh,p1);
suma(s,p1,s);
end;
writeln('Многочлен Ньютона ');
writeln(fi,'Многочлен Ньютона ');
vyvid(s);
fvyvid(s);
writeln;
writeln(fi);
w:=znach(s,t);
writeln('Значення в точцi iнтерполяцii=',w:5:3);
writeln(fi,'Значення в точцi iнтерполяцii=',w:5:3);
close(fi);
end.
Результат роботи програми
x y
3.00 4.00
7.00 10.00
11.00 22.00
15.00 26.00
19.00 33.00
Точка iнтерполяцii 13.000
Многочлен Ньютона
0.004x^ 4-0.183x^ 3+2.768x^ 2-14.087x+25.958
Значення в точцi iнтерполяцii=24.898
Виконання інтерполяції за допомогою пакету MathCad
ЧИСЕЛЬНЕ ІНТЕГРУВАННЯ ФУНКЦІЙ
Якщо функція f неперервна на відрізку [a;b] і відома її первісна F (F'(x)=f(x)), то справедлива формула Ньютона-Лейбніца:
.
Проте цією формулою важко і навіть практично неможливо скористатися тоді, коли первісну F не можна виразити в елементарних функціях або коли підінтегральну функцію задано графічним чи табличним способом.
У цих випадках треба будувати формули для наближеного обчислення визначених інтегралів. Особливо важливе значення мають методи чисельного інтегрування функцій, в яких для знаходження наближеного значення визначеного інтегралу використовуються значення підінтегральної функції та її похідних у скінченній кількості точок, що належать переважно проміжку інтегрування. Такі формули обчислення наближеного значення визначених інтегралів називають формулами механічних квадратур або квадратурними формулами.
В даному розділі розглянуті три методи чисельного інтегрування: метод прямокутників, метод трапецій, метод Сімпсона. Обчислено один визначений інтеграл всіма трьома методами. Для порівняння, те ж завдання виконане засобами пакету MathCad.
Метод прямокутників
Нехай потрібно обчислити інтеграл . Розіб'ємо ділянку інтегрування [a;b] на n рівних частин і помістимо точки, значення функції в яких входять до інтегральної суми, в лівих кінцях одержаних ділянок. Якщо вважати, що n достатньо велике, тобто довжина ділянок розбиття h= достатньо мала, то інтегральна сума повинна вже мало відрізнятися від величини інтегралу. Таким чином, ми одержимо наближену рівність:
h(
Ця формула називається формулою прямокутників (точніше, формулою лівих прямокутників). Якщо помістити точки, значення функції в яких входять до інтегральної суми, в правих кінцях одержаних ділянок, одержимо формулу правих прямокутників
h(
Якщо помістити точки, значення функції в яких входять до інтегральної суми, в серединах одержаних ділянок, одержимо формулу середніх прямокутників
h(
Програмна реалізація всіх трьох модифікацій методу практично однакова.
Завдання. Методом прямокутників знайти
Програмна реалізація
Здійснена на мові С.
Текст програми
/*програма знаходження iнтеграла методом лівих прямокутникiв*/
#include "stdio.h"
#include "math.h"
FILE *stream;
float f(float x)
{return(exp(1/x+sin(x)));};
main ()
{
/*опис змiнних*/
float a,b,h,s=0,x;
int n,i;
stream=fopen("respr.txt","w");
printf("Введiть a ");
scanf("%f",&a);
printf("введiть b ");
scanf("%f",&b);
fprintf(stream,"Iнтегрування методом прямокутникiв\n");
fprintf(stream,"Функцiя f(x)=exp(1/x+sin(x))\n");
fprintf(stream,"Межi a=%4.4f b=%4.4f\n",a,b);
n=10000;
h=(b-a)/n;
/*цикл по вузлах*/
for (i=0;i<=n;i=i++)
{
x=a+i*h;
s=s+f(x);
};
s=s*h;
/*вивiд результату*/
fprintf(stream,"Визначений iнтеграл=%4.3f\n",s);
fclose(stream);
}
Результат роботи програми
Iнтегрування методом прямокутникiв
Функцiя f(x)=exp(1/x+sin(x))
Межi a=5.0000 b=15.0000
Визначений iнтеграл=15.473
Метод трапецій
Замінивши прямокутники розбиття з попереднього розділу трапеціями, одержимо формулу трапецій:
h(
Завдання. Методом трапецій знайти
Програмна реалізація
Здійснена на мові С.
Текст програми
/*програма знаходження iнтеграла методом трапецiй*/
#include "stdio.h"
#include "math.h"
FILE *stream;
float f(float x)
{return(exp(1/x+sin(x)));};
main ()
{
/*опис змiнних*/
float a,b,h,s=0,x;
int n,i;
stream=fopen("restr.txt","w");
printf("Введiть a ");
scanf("%f",&a);
printf("введiть b ");
scanf("%f",&b);
fprintf(stream,"Iнтегрування методом трапецiй\n");
fprintf(stream,"Функцiя f(x)=exp(1/x+sin(x))\n");
fprintf(stream,"Межi a=%4.4f b=%4.4f\n",a,b);
n=10000;
h=(b-a)/n;
s=(f(a)+f(b))/2;
/*цикл по вузлах*/
for (i=1;i<=n-1;i=i++)
{
x=a+i*h;
s=s+f(x);
};
s=s*h;
/*вивiд результату*/
fprintf(stream,"s=%4.3f\n",s);
fclose(stream);
}
Результат роботи програми
Iнтегрування методом трапецiй
Функцiя f(x)=exp(1/x+sin(x))
Межi a=5.0000 b=15.0000
s=15.473
Метод Сімпсона
Замінивши прямолінійні трапеції з попереднього розділу криволінійними трапеціями, обмеженими параболами, і првівши відповідне інтегрування, одержимо формулу Сімпсона:
.
Формула Сімпсона більш точна, ніж формули прямокутників та трапецій. Це означає, що для досягнення тієї ж точності можна брати менше числоn ділянок розбиття, а при тому ж кроці h вона дає меншу абсолютну і відносну похибки.
Завдання. Методом Сімпсона знайти
Програмна реалізація
Здійснена на мові С.
Текст програми
/*програма знаходження iнтеграла методом Сiмпсона*/
#include "stdio.h"
#include "math.h"
float f(float x)
{return(exp(1/x+sin(x)));
};
main ()
{
/*опис змiнних*/
FILE *res;
float a,b,h,s=0,x;
int n,i;
res=fopen("ress.txt","w");
printf("введiть нижню межу iнтегрування");
scanf("%f",&a);
printf("введiть верхню межу iнтегрування");
scanf("%f",&b);
fprintf(res,"Iнтегрування ф-ii exp(1/x+sin(x)))\n");
fprintf(res,"Методом Сiмпсона\n");
fprintf(res,"a=%4.4f b=%4.4f\n",a,b);
n=10000;
h=(b-a)/n;
s=f(a)+f(b);
/*цикл по непарних вузлах*/
for (i=1;i<=n-1;i=i+2)
{
x=a+i*h;
s=s+4*f(x);
};
/*цикл по парних вузлах*/
for (i=2;i<=n-2;i=i+2)
{
x=a+i*h;
s=s+2*f(x);
};
s=s*h/3;
/*вивiд результату*/
fprintf(res,"s=%4.3f\n",s);
fclose(res);
}
Результат роботи програми
Iнтегрування ф-ii exp(1/x+sin(x)))
Методом Сiмпсона
a=5.0000 b=15.0000
s=15.473
Знаходження визначеного інтеграла засобами пакету MathCad
РОЗВ'ЯЗУВАННЯ ЗВИЧАЙНИХ ДИФЕРЕНЦІАЛЬНИХ РІВНЯНЬ
Часто задачі техніки і природознавства математично зводяться до відшукування розв'язку певного диференціального рівняння (або системи таких рівнянь), який задовольняє певні початкові умови (задача Коші). Далеко не завжди вдається проінтегрувати рівняння в скінченному вигляді. На практиці застосовують здебільшого наближене інтегрування диференціальних рівнянь. В даному розділі буде розглянуто три наближені методи розв'язку дифрівнянь і проілюстровано їх роботу на прикладі виконання завдання:
Знайти розв'язок диференціального рівняння
y'=y-x
який задовольняє початкову умову
y(0)=1.5
Це рівняння має загальний інтеграл
y=cex+x+1
Частковий розв'язок, який задовольняє початкову умову:
y= 0.5*ex+x+1
Метод Ейлера
Метод Ейлера застосовується для розв'язування диференціальних рівнянь виду y'=f(x,y), якщо відомі початкові умови y(x0)=y0 і крок h. Розв'язки диференціального рівняння знаходимо за ітераційною формулою:
yk+1=yk+f(xk,yk)*h
Програмна реалізація методу виконана на мові С.
Текст програми
/*програма числового розвязку дифрiвняння методом Ейлера*/
# include <stdio.h>;
# include <math.h>;
# include <conio.h>;
float fun(float a,float b) {
float c;
c=b-a;
return(c);
}
main() {
float k=1.5,h=0.25,x=0,y,yk;
FILE *res;
res=fopen("ejl.txt","w");
/*присвоення початкових значень*/
y=1.5;
/*yk-точне значення розвязку*/
/*початок циклу табуляцii наближеного i точного розвязку*/
while(x<=k) {
yk=0.5*exp(x)+x+1;
printf("x=%.2f y=%.2f точ.розв.=%.2f \n",x,y,yk);
fprintf(res,"x=%.2f y=%.2f точ.розв.=%.2f \n",x,y,yk);
y+=fun(x,y)*h;
x+=h;
};
fclose(res);
}
Результат роботи програми
x=0.00 y=1.50 точ.розв.=1.50
x=0.25 y=1.88 точ.розв.=1.89
x=0.50 y=2.28 точ.розв.=2.32
x=0.75 y=2.73 точ.розв.=2.81
x=1.00 y=3.22 точ.розв.=3.36
x=1.25 y=3.78 точ.розв.=4.00
x=1.50 y=4.41 точ.розв.=4.74
Уточнений метод Ейлера
Уточнений метод Ейлера застосовується для розв'язування диференціальних рівнянь y'=f(x,y), якщо відомі початкові умови y(x0)=y0 і крок h. Розв'язки диференціального рівняння знаходимо за ітераційною формулою:
yk+1=yk+f(xk+h/2,yk+f(xk,yk)*h/2)
Програмна реалізація методу виконана на мові С.
Текст програми
/*програма числового розвязку дифрiвняння методом Ейлера уточненим */
# include <stdio.h>
# include <math.h>
# include <conio.h>
float fun(float a,float b) { float c;
c=b-a;
return(c);
}
main() {
float k=1.5,h=0.25,x=0,y,x1,y1,x2,y2,yk;
FILE *res;
clrscr();
res=fopen("ejler.txt","w");
/*присвоення початкових значень*/
x1=0; y1=1.5;
/*yk-точне значення розвязку*/
yk=0.5*exp(x1)+x1+1;
printf("x Ейлер точний\n");
printf("x=%4.2f y=%4.4f yk=%4.4f\n",x1,y1,yk);
fprintf(res,"x Ейлер точний\n");
fprintf(res,"x=%4.2f y=%4.4f yk=%4.4f\n",x1,y1,yk);
/*початок циклу табуляцii наближеного i точного розвязку*/
while(x<k) {
/*подвiйний прорахунок*/
x2=x1+h/2;
x=x2+h/2;
y2=y1+fun(x1,y1)*h/2;
y=y2+fun(x2,y2)*h/2;
yk=0.5*exp(x)+x+1;
/*табуляцiя*/
printf("x=%4.2f y=%4.4f yk=%4.4f\n",x,y,yk);
fprintf(res,"x=%4.2f y=%4.4f yk=%4.4f\n",x,y,yk);
x1=x;
y1=y;
};
fclose(res);
}
Результат роботи програми
x Ейлер точний
x=0.00 y=1.5000 yk=1.5000
x=0.25 y=1.8828 yk=1.8920
x=0.50 y=2.3009 yk=2.3244
x=0.75 y=2.7636 yk=2.8085
x=1.00 y=3.2829 yk=3.3591
x=1.25 y=3.8737 yk=3.9952
x=1.50 y=4.5549 yk=4.7408
Метод Рунге-Кутта 4-го порядку
Цей метод служить для розв'язування диференціальних рівнянь виду y'=f(x,y), якщо відомі початкові умови y(x0)=y0 і крок h. Розв'язки диференціального рівняння знаходимо за ітераційними формулами:
yk+1=yk+
K1=hf(xk,yk)
K2=hf(
K3=hf(
K4=hf(xk+h,yk+K3)
Програмна реалізація методу виконана на мові С.
Текст програми
/*програма числового розвязку дифрiвняння*/
/*методом Рунге-Кута 4-го порядку*/
# include <stdio.h>
# include <math.h>
# include <conio.h>
float fun(float a,float b) {
float c;
c=b-a;
return(c);
}
main(){
float k=1.5,h=0.25,x,yk,y,k1,k2,k3,k4;
FILE *res;
clrscr();
res=fopen("rk4.txt","w");
/*присвоення початкових значень*/
x=0;
y=1.5;
/*yk-точне значення розвязку*/
yk=0.5*exp(x)+x+1;
printf("x Р-К точний\n");
printf("x=%4.2f z=%4.4f yk=%4.4f\n",x,y,yk);
fprintf(res,"x Р-К точний\n");
fprintf(res,"x=%4.2f z=%4.4f yk=%4.4f\n",x,y,yk);
/*початок циклу табуляцii наближеного i точного розвязку*/
while(x<k) {
k1=h*fun(x,y);
x+=h/2;
k2=h*fun(x,y+k1/2);
k3=h*fun(x,y+k2/2);
x=x+h/2;
k4=h*fun(x,y+k3);
y+=(k1+2*k2+2*k3+k4)/6;
yk=0.5*exp(x)+x+1;
/*табуляцiя*/
printf("x=%4.2f y=%4.4f yk=%4.4f\n",x,y,yk);
fprintf(res,"x=%4.2f y=%4.4f yk=%4.4f\n",x,y,yk);
};
fclose(res);
}
Результат роботи програми
x Р-К точний
x=0.00 z=1.5000 yk=1.5000
x=0.25 y=1.8920 yk=1.8920
x=0.50 y=2.3243 yk=2.3244
x=0.75 y=2.8085 yk=2.8085
x=1.00 y=3.3591 yk=3.3591
x=1.25 y=3.9951 yk=3.9952
x=1.50 y=4.7408 yk=4.7408
Подобные документы
Класичні та сучасні наближені методи розв'язання диференціальних рівнянь та їх систем. Класифікація наближених методів розв'язування. Розв'язування трансцендентних, алгебраїчних і диференціальних рівнянь, методи чисельного інтегрування і диференціювання.
отчет по практике [143,9 K], добавлен 02.03.2010Аналіз найвідоміших методів розв’язування звичайних диференціальних рівнянь і їх систем, користуючись рекомендованою літературою. Розробка відповідної схеми алгоритму. Розв’язання системи звичайних диференціальних рівнянь в за допомогою MathCAD.
лабораторная работа [412,4 K], добавлен 21.10.2014Розгляд теоретичних основ рівнянь з параметрами. Основні види даних рівнянь. Аналітичний та графічний методи розв’язування задач із використанням формул, властивостей функцій. Ознайомлення із системою розв’язування задач з параметрами для 9 класу.
курсовая работа [605,9 K], добавлен 29.04.2014Методи скінченних різниць або методи сіток як чисельні методи розв'язку інтегро-диференціальних рівнянь алгебри диференціального та інтегрального числення. порядок розв’язання задачі Діріхле для рівняння Лапласа методом сіток у прямокутної області.
курсовая работа [236,5 K], добавлен 11.06.2015Чисельні методи розв’язання систем нелінійних рівнянь: лінійні і нелінійні рівняння, метод простих ітерацій, метод Ньютона. Практичне використання методів та особливості розв’язання систем нелінійних рівнянь у пакеті Mathcad, Excel та на мові С++.
курсовая работа [2,0 M], добавлен 30.11.2010Вивчення методів розв'язання лінійної крайової задачі комбінуванням двох задач Коші. Переваги та недоліки інших методів: прицілювання, колокацій, Гальоркіна, найменших квадратів та ін. Пошук єдиного розв'язку звичайного диференціального рівняння.
курсовая работа [419,2 K], добавлен 29.08.2010Системи лінійних алгебраїчних рівнянь, головні означення. Коротка характеристика головних особливостей матричного способу, методу Жордано-Гаусса. Формули Крамера, теорема Кронекера-Капеллі. Практичний приклад розв’язання однорідної системи рівнянь.
курсовая работа [690,9 K], добавлен 25.04.2013Основні етапи розв'язування алгебраїчних рівнянь: аналіз задачі, пошук плану розв'язування та його здійснення; перевірка та розгляд інших способів виконання. Раціоналізація розв'язування алгебраїчних рівнянь вищих степенів методом заміни змінних.
курсовая работа [229,8 K], добавлен 13.05.2013Розв’язання систем лінійних рівнянь методом Жордана-Гауса. Еквівалентні перетворення системи, їх виконання як елемент методів розв’язування системи рівнянь. Базисні та вільні змінні. Лінійна та фундаментальна комбінації розв’язків, таблиці коефіцієнтів.
контрольная работа [170,2 K], добавлен 16.05.2010Розгляд найбільш відомих скінченно-різнецевих методів рішення рівнянь руху з непереривною силою: чисельна ітерація рівнянь Ньютона; алгоритм Бімана і Шофілда; метод Рунге-Кутта; методи Адамса, Крилова, Чаплигіна. Програма Рунге-Кутта на мові С#.
курсовая работа [359,5 K], добавлен 27.01.2011
