>

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;

:=

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;

:=

f

2

−

+

−

(

)

+

x

1

2

x

2

2

2

2 x

1

2

2 x

1

x

2

5

9

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;

:=

f

1

− + −

x

1

2

x

2

2

x

3

2

1

>

f[2]:=(x[1]-2)^2+(x[2]-2)^2+(x[3]-2)^2-36;

:=

f

2

+

+

−

(

)

−

x

1

2

2

(

)

−

x

2

2

2

(

)

−

x

3

2

2

36

>

f[3]:=x[2]-exp(-x[1]*x[3]);

:=

f

3

−

x

2

e

(

)

−x

1

x

3

>

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;

:=

f

1

−

x

1

(

)

sinh x

1

x

2

1
2

:=

f

2

−

+

−

(

)

+

x

1

2

x

2

2

2

2 x

1

2

2 x

1

x

2

5

9

10

>

for i to n do
for j to n do
A[i,j]:=diff(f[i],x[j]):
end do:
end do:

>

A;

















+

(

)

sinh x

1

x

2

x

1

(

)

cosh x

1

x

2

x

2

x

1

2

(

)

cosh x

1

x

2

−

+

4 (

)

+

x

1

2

x

2

2

x

1

4 x

1

2 x

2

5

+

4 (

)

+

x

1

2

x

2

2

x

2

10 x

1

x

2

4

>

w:=abs(f[1])<eps and abs(f[2])<eps;

w

<

−

x

1

(

)

sinh x

1

x

2

1
2

0.1000000000 10

-7

and

:=

<

−

+

−

(

)

+

x

1

2

x

2

2

2

2 x

1

2

2 x

1

x

2

5

9

10

0.1000000000 10

-7

>

x:=Vector([1.,1.]);

:=

x











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

( )

( )

( )

(

)

(

)

k

k

k

⋅

= −

A x

h

f x

)

(

)

(

)

1

(

k

k

k

h

x

x

+

=

+

>
>

x:=Vector([1.,1.]);

:=

x











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

)

(

)

(

)

(

)

1

(

k

k

k

x

f

B

x

x

⋅

−

=

+

>

x:=Vector([1.,1.]);

:=

x











1.
1.

>

B:=A^(-1);

:=

B











0.453736644776852105

-0.0388973461080574582

-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]};

:=

sys

{

}

,

−

+

−

(

)

+

x

1

2

x

2

2

2

2 x

1

2

2 x

1

x

2

5

9

10

−

x

1

(

)

sinh x

1

x

2

1
2

>

niew:={x[1],x[2]};

:=

niew

{

}

,

x

1

x

2

>

fsolve(sys,niew); # Maple wybiera pierwiastek

{

}

,

=

x

1

-1.516486440

=

x

2

0.2136586364

>

fsolve(sys,{x[1]=0..1,x[2]=0..1}); # my wybieramy pierwiastek

{

}

,

=

x

1

0.7613707931

=

x

2

0.8101727211

>
>