> restart:
Rozwiązywanie układów równań nieliniowych
Metoda graficzna
> with(plots):
Warning, the name changecoords has been redefined
Układ dwóch równań
> f[1]:=x[1]*sinh(x[1]*x[2])-1/2;
1
f1 := x1 sinh( x1 x2) -
2
> f[2]:=(x[1]^2+x[2]^2)^2-2*(x[1]^2-x[1]*x[2]^5)-9/10;
2
9
2 2 2 5
f2 := (x1 + x2 ) - 2 x1 + 2 x1 x2 -
10
> implicitplot([f[1]=0,f[2]=0],x[1]=-2..2,x[2]=-2..2,color=[blue,r
ed],numpoints=2000);
Układ trzech równań
> f[1]:=x[1]^2-x[2]^2+x[3]^2-1;
2 2 2
f1 := x1 - x2 + x3 - 1
> f[2]:=(x[1]-2)^2+(x[2]-2)^2+(x[3]-2)^2-36;
2 2 2
f2 := (x1 - 2) + ( x2 - 2 ) + ( x3 - 2 ) - 36
> f[3]:=x[2]-exp(-x[1]*x[3]);
( -x x )
1 3
f3 := x2 - e
> implicitplot3d([f[1]=0,f[2]=0,f[3]=0],x[1]=-5..8,x[2]=-5..8,x[3]
=-5..8,color=[blue,green,red],orientation=[65,135],numpoints=500
0);
>
Metoda Newtona-Raphsona
>
Układ dwóch równań
> with(LinearAlgebra):
> n:=2;
n := 2
> eps:=10.^(2-Digits);
eps := 0.1000000000 10-7
> f:=Vector(n):A:=Matrix(n):
> f[1]:=x[1]*sinh(x[1]*x[2])-1/2;
f[2]:=(x[1]^2+x[2]^2)^2-2*(x[1]^2-x[1]*x[2]^5)-9/10;
1
f1 := x1 sinh( x1 x2) -
2
2
9
2 2 2 5
f2 := (x1 + x2 ) - 2 x1 + 2 x1 x2 -
10
> for i to n do
for j to n do
A[i,j]:=diff(f[i],x[j]):
end do:
end do:
> A;
2
îÅ‚sinh(x x2 ) + x1 cosh( x1 x2 ) x2 Å‚Å‚
x1 cosh( x1 x2)
ïÅ‚ śł
1
ïÅ‚ śł
2 2 5 2 2 4
ïÅ‚
ïÅ‚ śł
4( x1 + x2 ) x1 - 4 x1 + 2 x2 4( x1 + x2 ) x2 + 10 x1 x2 śł
ðÅ‚ ûÅ‚
> w:=abs(f[1])1
w := x1 sinh( x1 x2) - < 0.1000000000 10-7 and
2
2
9
2 2 2 5
( x1 + x2 ) - 2 x1 + 2 x1 x2 - < 0.1000000000 10-7
10
> x:=Vector([1.,1.]);
îÅ‚ 1. Å‚Å‚
ïÅ‚ śł
x := ïÅ‚ śł
ðÅ‚ 1. ûÅ‚
> i:=0;
i := 0
> while not w do
x:=x-A^(-1).f:
i:=i+1:
end do:
> x;i;
îÅ‚0.761370793108574806Å‚Å‚
ïÅ‚ śł
ïÅ‚ śł
ðÅ‚0.810172721149395758ûÅ‚
5
> f;
îÅ‚ 0. Å‚Å‚
ïÅ‚ śł
ïÅ‚ śł
ðÅ‚0.14 10-8ûÅ‚
>
Paktyczne implementacje metody Newtona-Raphsona
>
1. Liniowy ukad równań ze względu na elementy wektora h
A(x(k )) Å"h(k ) = -f (x(k ))
x(k +1) = x(k) + h(k)
>
> x:=Vector([1.,1.]);
îÅ‚ 1. Å‚Å‚
ïÅ‚ śł
x := ïÅ‚ śł
ðÅ‚ 1. ûÅ‚
> i:=0;
i := 0
> while not w do
h:=eval(LinearSolve(A,-f)):
x:=x+h:
i:=i+1:
end do:
> x;i;
îÅ‚0.7613707930Å‚Å‚
ïÅ‚ śł
ïÅ‚ śł
ðÅ‚0.8101727213ûÅ‚
5
> h;
îÅ‚
0.2377202314 10-6 Å‚Å‚
ïÅ‚ śł
ïÅ‚ śł
ðÅ‚-0.3330045624 10-6ûÅ‚
> f[1];f[2];
0.
0.22 10-8
>
1. Iteracja ze stałą macierzą B
x(k +1) = x(k) - B Å" f (x(k))
> x:=Vector([1.,1.]);
îÅ‚ 1. Å‚Å‚
ïÅ‚ śł
x := ïÅ‚ śł
ðÅ‚ 1. ûÅ‚
> B:=A^(-1);
îÅ‚ 0.453736644776852105 -0.0388973461080574582Å‚Å‚
ïÅ‚ śł
B := ïÅ‚ śł
ðÅ‚-0.151245548258950702 0.0685213375915747076 ûÅ‚
> i:=0;
i := 0
> while not w do
x:=x-B.f:
i:=i+1:
end do:
> x;i;
îÅ‚0.761370792799353380Å‚Å‚
ïÅ‚ śł
ïÅ‚ śł
ðÅ‚0.810172722281518598ûÅ‚
35
> f[1];f[2];
0.5 10-9
0.82 10-8
>
RozwiÄ…zanie Maple'a
> x:='x';
x := x
> sys:={f[1],f[2]};
2
9 1
2 2 2 5
sys := { ( x1 + x2 ) - 2 x1 + 2 x1 x2 - , x1 sinh( x1 x2) - }
10 2
> niew:={x[1],x[2]};
niew := {x1, x2 }
> fsolve(sys,niew); # Maple wybiera pierwiastek
{x1 = -1.516486440, x2 = 0.2136586364 }
> fsolve(sys,{x[1]=0..1,x[2]=0..1}); # my wybieramy pierwiastek
{ x1 = 0.7613707931, x2 = 0.8101727211 }
>
>