Laboratorium Metod Numerycznych

Ćw. 5

Temat: Numeryczne rozwiązywanie równań i układów równań nieliniowych

Ćwikliński Radosław 163509

Majer Łukasz 163557

1. Cel ćwiczenia

Praktyczne sprawdzenie wiedzy n/t popularnych metod iteracyjnych rozwiazywania równań i układów równań nieliniowych. Porównanie przydatności poszczególnych metod do wyznaczania zer funkcji określonych typów. Prześledzenie związku miedzy rzędem metody iteracyjnej a szybkością zbieżności ciągu przybliżeń.

2. Porównanie metody bisekcji, siecznych i Newtona- Raphsona

Funkcja 1: f(x)=arctg(4x)-πx

Pochodna funkcji 1: f'(x)=-π + 4/16x²+1

Wartość dokładna: x_* = 0.25

Punkt startowy dla metody Newtona: x₀ = 0.3

Punkt startowy dla metody siecznych x₀ = 0.3 x₁ = 0.29

Punkt startowy dla metody bisekcji x_a = 0.225 x_b = 0.3

M-pliki:

-metoda bisekcji

x=[-2:0.01:6];

y=[(x.^2)-2];

plot(x,y),grid;

b=(0.3);

a=(0.225);

xi=a;

for i=1:1:20

err(i,:)=[0.25-xi];

hold on

plot(a,0,'r*',b,0,'b^')

if (((atan(4*b)-pi*b)*(atan(4*xi)-pi*xi)) < 0)

a=xi;

xi=(b+a)/2

else

b=xi;

xi=(b+a)/2

end

end

figure

hold on

plot(xi,0,'b*'),grid

figure

zak_x=[1:20];

plot(zak_x,err)

-metoda siecznych

x(1)=0.3;

x(2)= 0.29;

c=0.25

for i=2:1:5

fxi=atan(4*x(i))-pi*x(i);

fxi_1=atan(4*x(i-1))-pi*x(i-1);

x(i+1)=x(i) - (fxi * (x(i) - x(i-1)))/(fxi - fxi_1);

blad(i)=c-x(i+1);

end

plot(x,'b*'),grid,xlabel('numer przyblizenia'),ylabel('wartosc');

plot(blad,'r-'),grid,xlabel('numer przyblizenia'),ylabel('blad');

- metoda Newtona- Raphsona

x(1)=0.3;

c=0.25

for i=1:1:4

fx=atan(4*x(i))-pi*x(i);

fxx=-pi+(4/(1+16*x(i)^2));

x(i+1)=([x(i)]-[fx/fxx])

blad(i)=c-x(i+1);

end

plot(x,'b*'),grid,xlabel('numer przyblizenia'),ylabel('wartosc')

plot(blad,'b*'),grid,xlabel('numer przyblizenia'),ylabel('blad')

Funkcja 2: f(x)=2(cos(πx)+ πx)- π

Pochodna funkcji 1: f'(x)= 2 (π-π sin(π x))

Wartość dokładna: x_* = 0.5

Punkt startowy dla metody Newtona: x₀ = 0.45

Punkt startowy dla metody siecznych x₀ = 0.45 x₁ = 0.46

Punkt startowy dla metody bisekcji x_a = 0.45 x_b = 0.6

M-pliki:

-metoda bisekcji

x=[0:0.01:4];

y=[2*(cos(pi*x)+pi*x)-pi];

plot(x,y),grid;

b=(0.45);

a=(0.6);

xi=a;

for i=1:1:20

err(i,:)=[0.5-xi];

hold on

plot(a,0,'r*',b,0,'b^')

if (((2*(cos(pi*b)+pi*b)-pi)*(2*(cos(pi*xi)+pi*xi)-pi)) < 0)

a=xi;

xi=(b+a)/2

else

b=xi;

xi=(b+a)/2

end

end

figure

hold on

plot(xi,0,'b*'),grid

figure

zak_x=[1:20];

plot(zak_x,err)

Wykres błędu:

-metoda siecznych

x(1)=0.45;

x(2)= 0.46;

c=0.25

for i=2:1:8

fxi=2*(cos(pi*x(i))+pi*x(i))-pi;

fxi_1=2*(cos(pi*x(i-1))+pi*x(i-1))-pi;

x(i+1)=x(i) - (fxi * (x(i) - x(i-1)))/(fxi - fxi_1);

blad(i)=c-x(i+1);

end

plot(blad,'b*'),grid,xlabel('numer przyblizenia'),ylabel('blad');

-metoda Newtona-Raphsona

x(1)=0.45;

c=0.5

for i=1:1:4

fx=2*(cos(pi*x(i))+pi*x(i))-pi;

fxx=2*(pi-pi*sin(pi*x(i)))

x(i+1)=([x(i)]-[fx/fxx])

blad(i)=c-x(i+1);

end

plot(blad,'b*'),grid,xlabel('numer przyblizenia'),ylabel('blad');

3. Porównanie pierwiastków wielomianu otrzymanych różnymi metodami:

P=[1 -8 12 50 -250 318 696 -2394 585 9050 -20500 15000];

Pierwiastki otrzymane za pomocą funkcji roots:

5.0000

-3.0000

-2.0000 + 1.0000i

-2.0000 - 1.0000i

1.0000 + 2.0000i

1.0000 - 2.0000i

1.0000 + 2.0000i

1.0000 - 2.0000i

2.0000 + 0.0000i

2.0000 - 0.0000i

2.0000

Pierwiastki otrzymane za pomocą funkcji laguerre.m:

1.9995

1.9995

2.0010

0.9999 - 2.0009i

-2.0000 - 1.0000i

1.0003 + 2.0007i

-2.0000 + 1.0000i

1.0001 - 1.9991i

0.9997 + 1.9993i

-3.0000 - 0.0000i

5.0000 - 0.0000i

Pierwiastki otrzymane za pomocą funkcji bairstow.m:

1.9996 - 0.0004i

1.9996 + 0.0004i

-3.0003

2.0007

Pierwiastki otrzymane za pomocą funkcji lehmer.m:

1.9999 + 0.0000i

4.9999 + 0.0000i

2.1217 + 0.0000i

1.8589 + 0.0000i

0.9655 + 2.0582i

1.0431 + 1.9385i

-1.9972 + 1.0006i

-3.0004 + 0.0013i

-2.0016 - 1.0004i

0.9489 - 2.0761i

1.0611 - 1.9220i

4. Rozwiązanie układu równań za pomocą metody Newtona

Funkcje:

Punkt startowy: (x₀,y₀) = (3.3; 0.8)

Wartość poprawna (x_*,y_*) = (3 1)

function funkcje

eps=1*10^-9;

syms x y;

f1=inline('x*y-x+y-1','x','y');

f2=inline('x^2+2*y^2-11','x','y');

F=[f1(x,y); f2(x,y)];

J=jacobian(F,[x y]);

F=inline(F);

J=inline(J);

X=[3.3 0.8 ]';

F(X(1),X(2));

J(X(1),X(2));

X1=X-inv(J(X(1),X(2)))*F(X(1),X(2));

for k=1:100

X1=X-J(X(1),X(2))\F(X(1),X(2));

if(norm(X1-X)<eps)

break

end

X=X1;

end

punkt=X1

punkt =

3

1

---