> restart:
Metody przybliżone rozwiązywania równań nieliniowych Znaleźć pierwiastek równania w przedziale [0, 1]: x
e 3 x 0
>
Metoda połowienia przedziału (bisekcji)
>
> f:=x->exp(x)-3*x;
f := xe
x
3 x
> plot(f(x),x=0..1);
> a:=0.;b:=1.;eps:=1.*10^(1-Digits);
:=
a
0.
:=
b
1.
eps := 0.1000000000 10-8
> i:=0;
:=
i
0
> while b-a>eps do
x:=(a+b)/2:
if f(x)=0 then break elif f(a)*f(x)<0
then b:=x
else a:=x
end if;
i:=i+1:
end do:
> xbis:=x; bis:=i; # wartość pierwiastka i liczba iteracji
:=
xbis
0.6190612875
:=
bis
27
> x:='x':
Metoda iteracji prostej
> a:=0.;b:=1.;
:=
a
0.
:=
b
1.
> g:=x->1/3*exp(x); # g:=x->ln(3*x); 1
:=
g
x
e x
3
> gp:=D(g);
1
:=
gp
x e x
3
> plot([abs(gp(x)),1],x=a..b); # plot([abs(gp(x)),1],x=b..2*b);
> x:=a; # x:=b;
x := 0.
> i:=0;
:=
i
0
x g( x)
> while f(x)*f(x-eps)>0 and f(x)*f(x+eps)>0 do x:=g(x);
i:=i+1:
end do:
> xip:=x;ip:=i;
xip := 0.6190612853
:=
ip
40
>
Metoda Newtona (stycznych)
>
> x:='x':
> a:=0.; b:=1.;
:=
a
0.
:=
b
1.
> fp:=D(f);
fp := xe
x
3
> fb:=D(fp);
fb := exp
> if f(a)*fb(a)>0 then xold:=a else xold:=b end if; xold := 0.
f ( x )
k
x
x
k 1
k
f (
x )
k
> xnew:=xold-f(xold)/fp(xold); xnew := 0.5000000000
> i:=0;
:=
i
0
f ( x )
k
x
x
k 1
k
f (
x )
k
> while abs(xnew-xold)>eps do xold:=xnew:
xnew:=xold-f(xold)/fp(xold); i:=i+1:
> xNewt:=xnew; Newt:=i; # wartość pierwiastka i liczba iteracji xNewt := 0.6190612866
Newt := 5
>
Metoda Halley'a
>
> x:='x':
> if f(a)*fb(a)>0 then xold:=a else xold:=b end if; xold := 0.
> xnew:=xold-2*f(xold)*fp(xold)/(2*fp(xold)^2-fb(xold)*f(xold)); xnew := 0.5714285714
> i:=0;
:=
i
0
2 f ( x ) f (
x )
k
k
x
x
k 1
k
2
2 f (
x ) f ( x ) f ( x ) k
k
k
> while abs(xnew-xold)>eps do xold:=xnew:
xnew:=xold-2*f(xold)*fp(xold)/(2*fp(xold)^2-fb(xold)*f(xold)); i:=i+1:
end do:
> xHall:=xnew; Hall:=i; # wartość pierwiastka i liczba iteracji xHall := 0.6190612869
Hall := 3
>
Porównanie dokładności i liczby iteracji wszystkich metod
> 'xbis'=xbis,bis;'xip'=xip,ip;'xNewt'=xNewt,Newt;'xHall'=xHall,Hall;
xbis
0.6190612875, 27
xip0.6190612853, 40
xNewt0.6190612866, 5
xHall0.6190612869, 3
>
Komenda w Maple'u
> x:='x':
> f(x);
e
x
3 x
> xMaple:=fsolve(f(x)=0,x=0..1); xMaple := 0.6190612867
>
Zmiana dokładności rozwiązania
> Digits:=50; x:='x':eps:=1.*10^(1-Digits);
>
Zmiana przedziału [1, 2]