Реферат: Обчислювальні методи розв’язку нелінійних рівнянь

ЗМІСТ

Обчислювальні методи розв’язку нелінійних рівнянь 3

Методи хорд, дотичних, хорд-дотичних 3

Метод хорд 6

Метод дотичних 7

Метод хорд-дотичних 8

Метод ітерацій 8

Завдання. 9

Програмна реалізація методу 9

Метод половинного поділу 10

Завдання. 10

Програмна реалізація методу 10

МЕТОДИ ЛІНІЙНОЇ АЛГЕБРИ 12

Знаходження визначника матриці методом алгебраїчних доповнень 12

Розв’язок системи рівнянь методом Крамера 14

Розв’язок системи рівнянь методом Гауса 17

20

Знаходження визначника матриці методом Гауса 20

Знаходження оберненої матриці за допомогою 22

Розв’язок системи лінійних рівнянь 25

за допомогою методу Зейделя 25

^ ІНТЕРПОЛЮВАННЯ ФУНКЦІЙ 27

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

Інтерполяційний многочлен Ньютона 30

ЧИСЕЛЬНЕ ІНТЕГРУВАННЯ ФУНКЦІЙ 35

Метод прямокутників 35

Метод трапецій 37

Метод Сімпсона 37

РОЗВ’ЯЗУВАННЯ ЗВИЧАЙНИХ ДИФЕРЕНЦІАЛЬНИХ РІВНЯНЬ 40

Метод Ейлера 40

Уточнений метод Ейлера 41

Метод Рунге-Кутта 4-го порядку 42
^ Обчислювальні методи розв’язку нелінійних рівнянь

Нехай задано рівняння з однією змінною

f(x)=0, (1)

де функція f(x) визначена і неперервна на деякому проміжку .

^ Розв’язати рівняння - означає знайти множину його коренів, тобто таких значень x, при яких рівняння (1) перетворюється в тотожність. Корінь рівняння (1) називають ще нулем функції f(x). Якщо функція f(x) - алгебраїчний многочлен, то рівняння (1) називають алгебраїчним. Якщо функція f(x) містить тригонометричні, показникові, логарифмічні і т.д. функції, то рівняння (1) називають трансцендентним (нелінійним).

Знайти точні значення коренів заданого рівняння можна лише для найпростіших функцій f(x); алгебраїчних многочленiв не вище четвертого степеня, деяких многочленiв степеня n5 і деяких трансцендентних функцій.

Універсальних методів для знаходження точних значень коренів алгебраїчних рівнянь степеня n5 і трансцендентних рівнянь не існує. Тому важливого значення набувають наближені методи знаходження коренів з достатньою для практики точністю. Задача знаходження коренів рівняння (1) вважається розв’язаною, якщо корені обчислені із наперед заданою точністю.

Знаходження наближених коренів рівняння (1) складається з двох етапів:

1) відокремлення коренів, тобто знаходження досить малих відрізків, на кожному з яких міститься один і тільки один корінь рівняння;

2) обчислення коренів з наперед заданою точністю.

До всіх представлених методів розроблені програмні реалізації на мові С. Для порівняння, завдання виконані також за допомогою пакету MathCad.
^ Методи хорд, дотичних, хорд-дотичних
В даному розділі реалізовані такі методи
I.метод хорд II.метод дотичних III.метод хорд-дотичних
Ці методи полягають у наближеній заміні на досить малому відрізку функції f(x) лінійною функцією, яка відповідає хорді, проведеній через дві задані точки (метод хорд), дотичною, проведеною в заданій точці (метод дотичних), або ж ці методи комбінуються (метод хорд-дотичних).


Метод хорд

Якщо на інтервалі [a;b] неперервна функція F(x) задовольняє умову F(a)*F(b)<0, то корінь рівняння F(x)=0 наближено знаходиться за рекурентною формулою

xn=xn-1- (2)

Похибка обчислюється за формулою

|xn-xn-1|< (3)


Г
Xn-1

Xn

Xn-2
рафічна ілюстрація приведена на рис.1.


^ Метод дотичних
Корінь рівняння F(x)=0 обчислюється за ітераційною формулою

xk+1=xk- (4)


Графічна ілюстрація приведена на рис.2.










Xk+1






Xk

рис.2

Метод хорд-дотичних

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






an+1







bn



bn+1



an







рис.3

Завдання.

Використовуючи методи хорд, дотичних, хорд-дотичних знайти корінь функції f(x)=x3-2x2+x-3 на інтервалі [2,1;2,2].

Програмна реалізація методів
Метод хорд

Текст програми


#include

#include

/*метод хорд*/


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

#include

/*метод дотичних*/


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

#include

/*метод хорд-дотичних*/


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).

Нехай виділений інтервал ізоляції кореня цього рівняння, і x0 - довільна точка цього інтервалу (нульове наближення). Для одержання наступного наближення в праву частину рівняння (4) замість x підставляємо x0 і т.д. Ітераційний процес збігається, якщо відображення (x) - стягуюче.
Завдання.
За допомогою методу ітерацій знайти корінь функції f(x)= при нульовому наближенні 2,1.
Програмна реалізація методу

Текст програми


#include

#include

/*метод 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) задовольняє такі умови: неперервність і диференційованість на відрізку , набуває на кінцях цього відрізка значень різних знаків, а похідна f’(x) зберігає сталий знак всередині відрізка . Уточнення інтервалу відбувається методом половинного поділу.
Завдання.
За допомогою методу половинного поділу знайти корінь функції f(x)=x3-x2+1 на відрізку [2;3].
^ Програмна реалізація методу

Текст програми

/*пошук коpенiв методом половинного подiлу*/

#include

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

#include

/*п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
if (fmod(j,2)==0) zn=1;

else zn=-1;

k=1;l=0;n1=n-1;

for (k1=0;k1
l=0;

for (l1=0;l1
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

for(j=0;j

scanf("%f",&a);

mas[i][j]=a;

};

};

fprintf(res,"Матриця\n");

for(i=0;i


for(j=0;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

#include

/*п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
if (fmod(j,2)==0) zn=1;

else zn=-1;

k=1;l=0;n1=n-1;

for (k1=0;k1
l=0;

for (l1=0;l1
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

for (j=0;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


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

for (j=0;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

for (k=0;k

for(j=0;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

s=0;

for (j=0;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)
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

#include

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

#include

/* 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
printf("Введiть значення y ");

for(i=0;i
printf("Введiть значення т-ки iнтерполяцii");

scanf("%f",&t);

for (i=0;i
l=0;

for(i=0;i
{

p=1;

for(j=0;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

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
{puts("ввiд x");

scanf("%f",&x[i]);

fprintf(stream,"%5.2f ",x[i]);

};

fprintf(stream,"\ny ");

for (i=0;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
for(i=0;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
{ 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 вона дає меншу абсолютну і відносну похибки.

З