Rješavanje sustava nelinearnih jednadžbi

Primjer 1

f1[x1_,x2_]:=x2  (x1-1)-1
f2[x1_,x2_]:=x1^2-x2^2-1
ContourPlot[{f1[x1,x2]==0, f2[x1,x2]==0},{x1,-3,3},{x2,-3,3},PlotRange->{-2,2}, AspectRatio->.7]

Primjer 2

f1[x1_,x2_]:=x1^2+x2^2-4
f2[x1_,x2_]:=x1^2+9x2^2-9
ContourPlot[{f1[x1,x2]==0, f2[x1,x2]==0},{x1,-3,3},{x2,-3,3},PlotRange->{-2,2}, AspectRatio->.7]

Primjer 3

f1[x1_,x2_]:=x1^2-x2+al
f2[x1_,x2_]:=-x1+x2^2+al

al=1;
sl1=ContourPlot[{f1[x1,x2]==0, f2[x1,x2]==0},{x1,-3,3},{x2,-3,3},PlotRange->{-2,2}, AspectRatio->.7 ,PlotLabel->"al=1"];
al=1/4;
sl2=ContourPlot[{f1[x1,x2]==0, f2[x1,x2]==0},{x1,-3,3},{x2,-3,3},PlotRange->{-2,2}, AspectRatio->.7 ,PlotLabel->"al=1/4"];
al=0;
sl3=ContourPlot[{f1[x1,x2]==0, f2[x1,x2]==0},{x1,-3,3},{x2,-3,3},PlotRange->{-2,2}, AspectRatio->.7,PlotLabel->"al=0" ];
al=-1;
sl4=ContourPlot[{f1[x1,x2]==0, f2[x1,x2]==0},{x1,-3,3},{x2,-3,3},PlotRange->{-2,2}, AspectRatio->.7 ,PlotLabel->"al=-1"];
Show[GraphicsArray[{{sl1,sl2},{sl3,sl4}}]]

Newtonova metoda

Definiranje funkcije i Jacobijana

Clear[f,J,x1,x2];n=2;
pom[x_]:={f1[x[1],x[2]],f2[x[1],x[2]]};
f[x]:=pom[x];
J[x_]:=Table[D[f[x][[i]],x[j]],{i,n},{j,n}]
Print["f(x)=",MatrixForm[f[x]/.{x[1]->x1,x[2]->x2}],",   J(x)=",MatrixForm[J[x]/.{x[1]->x1,x[2]->x2}]]
ContourPlot[{f1[x1,x2]==0, f2[x1,x2]==0},{x1,-3,3},{x2,-3,3},PlotRange->{-2,2}, AspectRatio->.7]

Pocetna aproksimacija

Iterativni proces

w=Table[x[i]->x0[[i]],{i,n}];
a=J[x]/. w ;
b=-f[x]/.w;
Print["x0=",N[x0],"  f0= ",-N[b]]
s0=LinearSolve[a,b]; x1=x0+s0; x0=x1;

Modul za Newtonovu metodu

SNJ[n_,x0_,eps_,it_,f_]:=Module[{a,b,s0,xs,xn,k=0},w=Table[x[i]->x0[[i]],{i,n}];
J[x_]:=Table[D[f[x][[i]],x[j]],{i,n},{j,n}];
a=J[x]/. w ;
b=-f[x]/.w;
xs=x0;
Print["x0=",N[x0],"   f0=",-N[b]];
             While[k=k+1;Max[Abs[b]]>eps&&k<it, 
                          s0=LinearSolve[a,b];
                          xn=xs+s0; 
                          xs=xn; 
                          w=Table[x[i]->xs[[i]],{i,n}];
                          a=J[x]/. w ;
                          b=-f[x]/.w;
                          Print["x",k,"=",N[xn],"   f=",-N[b]];
                        ]
         ]

Poziv modula "Newtonova metoda"

SNJ[2,{2,0},.00005,10,f]

Quasi-Newtonova metoda (Broydenova metoda- 1965)

Ulazni podaci

B0=H0=IdentityMatrix[n]; x0={2.,2.};k=0;it=100;eps=.00005;
Print["x0=",x0,"   f0=",f[x]/.Table[x[i]->x0[[i]],{i,n}]]

Broydenova metoda

While[k=k+1;
w0=Table[x[i]->x0[[i]],{i,n}];
s0=-Inverse[B0].f[x]/.w0;
xn=x0+s0;
w1=Table[x[i]->xn[[i]],{i,n}];
fn=f[x]/.w1; f0=f[x]/.w0; y=fn-f0;
Max[Abs[fn]]>eps&&k<it, 
Print["x",k,"=",xn,"     f",k,"=",fn]
If[Apply[Plus,s0 Inverse[B0].y]==0,B1=B0, v=s0/Apply[Plus,s0  s0]];
vv=y-B0.s0;B1=B0+Table[vv[[i]] v[[j]],{i,n},{j,n}];
x0=xn; B0=B1]

Inverzna Broydenova metoda

While[k=k+1;
w0=Table[x[i]->x0[[i]],{i,n}];
s0=-H0.f[x]/.w0;
xn=x0+s0;
w1=Table[x[i]->xn[[i]],{i,n}];
fn=f[x]/.w1; f0=f[x]/.w0; y=fn-f0;
Max[Abs[fn]]>eps&&k<it, 
Print["x",k,"=",xn,"     f",k,"=",fn]
If[y==Table[0,{i,n}],H1=H0, v=(Transpose[H0].s0)/Apply[Plus,(H0.y)  s0]];
vv=s0-H0.y;H1=H0+Table[vv[[i]] v[[j]],{i,n},{j,n}];
x0=xn; H0=H1]