solving non linear equations

clc
clear all 
syms  x y z xn xnp
double err 

ans = 1×3

101 114 114

int16 n;
err=10^-4

err = 1.0000e-04

n=2;
f1= [ x^2+x-y^2+x*y+z^2-3
 5*(y^2)+y+y*z-x^2+z-7
  x+y+z+y*z-z^3-3]

f1 = 

g=jacobian([f1],[x,y,z])

g = 

t=inv(g)

t = 

xn=sym([0.5;0.5;0.5]) 

xn = 

xnp=-xn

xnp = 

i=0

i = 0

while max(abs(xnp-xn))> err 
    xn=xnp
   fc=f1
   jc=t
   fc(xn)=subs(fc,[x;y;z],xn)
   jc(xn)=subs(jc,[x;y;z],xn)
   xnp=xn-(fc*jc)
   fc=[]
   jc=[]
    
    fprintf('Iteration %d: x=%.18f',i, xnp); 
   i=i+1 
end

xn = 

fc = 

jc = 

Array indices must be positive integers or logical values.

Error in sym/privsubsasgn (line 1311)
L_tilde2 = builtin('subsasgn',L_tilde,struct('type','()','subs',{varargin}),R_tilde);

Error in indexing (line 1142)
C = privsubsasgn(L,R,inds{:});

i am tryin here to write my own code to solve non linear system based on Newton method am pretty sure about the algorithm and how to use the mehtod in solving non linear system however am not sure about the syntax of matlab coding . can anyone please help with an explanation of these errors and how to avoid them. thanks in advance .

Answer 1

[VBBV]

fc(:,i+1)=subs(fc,[x;y;z],xn)
 jc(:,i+1)=subs(jc,[x;y;z],xn)

[VBBV]

Use the increment i for arrays. Matlab uses 1 based indexing for storing elements in array. In your case xn is a array with non integers or fractions which you tried to use as index to store the result in fc and jc variables

[VBBV]

clc
clear all 
syms  x y z xn xnp
double err; 
int16 n;
err=10;
n=2;
f1= [ x^2+x-y^2+x*y+z^2-3
 5*(y^2)+y+y*z-x^2+z-7
  x+y+z+y*z-z^3-3];

g=jacobian([f1],[x,y,z]);
t=inv(g);
xn=sym([0.5;0.5;0.5]) ;
xnp=-xn;
i=0;
 
while max(abs(xnp-xn)) < err      % check the condition 
    xn=xnp;
   fc=f1;
   jc=t;
   fc(:,i+1)=subs(fc,[x;y;z],xn);
   jc(:,:,i+1)=subs(jc,[x;y;z],xn); %  use the iteration index i
   xnp=xn-(fc(:,i+1)'*jc(:,:,i+1)).';
 
    fprintf('Iteration %d:\n',i+1); 
    fprintf(' x=%.18f\n',xnp)
   i=i+1; 
end

Iteration 1:

 x=-6.092592592592592560
 x=-2.388888888888888840
 x=1.592592592592592560

Iteration 2:

 x=-4.350537376611543827
 x=-3.586670317708851030
 x=0.423387958776677398

Iteration 3:

 x=-2.715600164121943827
 x=-2.647869221303980947
 x=-3.786070857048641436

Iteration 4:

 x=0.460412999441244641
 x=-1.594645782058791594
 x=-3.195679571515712247

Iteration 5:

 x=59.497556895760055795
 x=11.977820988089051824
 x=-14.562217932844697899

[amr makhlouf]

may i ask why did you make the error 10 instead of 10e-4 ? and use the < err isntead of > err ?

Answer 2

[Torsten]

syms  x y z  
errX = 10;
errF = 10;
imax = 25;
TolX = 1e-8;
TolF = 1e-8;
f = [ x^2+x-y^2+x*y+z^2-3
      5*(y^2)+y+y*z-x^2+z-7
      x+y+z+y*z-z^3-3];

J = jacobian(f,[x,y,z]);
Jinv = inv(J);
xi = [0.5;0.5;0.5] ;
i = 0;
 
while (errX > TolX  || errF > TolF) &&   i < imax    % check the condition 
    fi = double(subs(f,[x y z],[xi(1) xi(2) xi(3)]));
    Ji = double(subs(Jinv,[x y z],[xi(1) xi(2) xi(3)]));
    xip1 = xi - Ji*fi;
    i = i + 1;
    errX = max(abs(xip1-xi))
    errF = norm(fi)
    fprintf('Iteration %d:\n',i); 
    fprintf(' x=%.18f\n',xi)
    fprintf(' f=%.18f\n',fi)
    xi = xip1;
end

errX = 5.6707

errF = 5.3341

Iteration 1:

 x=0.500000000000000000
 x=0.500000000000000000
 x=0.500000000000000000

 f=-2.000000000000000000
 f=-4.750000000000000000
 f=-1.375000000000000000

errX = 6.2083

errF = 237.9880

Iteration 2:

 x=-1.134146341463414753
 x=-0.329268292682926678
 x=6.170731707317072434

 f=35.495092207019631303
 f=-3.934562760261748782
 f=-235.293198009315005947

errX = 2.8092

errF = 75.8640

Iteration 3:

 x=5.074114328386549744
 x=-2.572548251018186072
 x=4.011835134589464502

 f=24.244163248276837663
 f=-8.537966273744203249
 f=-71.377007209879081984

errX = 1.0865

errF = 22.7986

Iteration 4:

 x=2.264894613092790721
 x=-1.503611390419629990
 x=2.631309321385490119

 f=5.652062414573667759
 f=-3.654280277805187271
 f=-21.782503986816859509

errX = 0.8495

errF = 7.6660

Iteration 5:

 x=3.042354881676863965
 x=-2.084842561940826400
 x=1.544790017024043172

 f=0.995254850508691113
 f=1.716222792667095653
 f=-7.404791765514490365

errX = 1.0762

errF = 2.8990

Iteration 6:

 x=3.269055297861723108
 x=-2.016092559757794778
 x=0.695304144141349112

 f=0.783878418233097340
 f=-0.086162420312693569
 f=-2.789673925096745943

errX = 1.2439

errF = 1.6680

Iteration 7:

 x=3.376106628947205035
 x=-1.999442692755602824
 x=-0.380902932842055564

 f=1.171186832390128885
 f=-0.027992601829813279
 f=-1.187381330258901135

errX = 0.3648

errF = 3.7617

Iteration 8:

 x=2.565025233115954428
 x=-1.654064465972460596
 x=-1.624811290170253830

 f=1.805745058289646821
 f=-0.491041295718978865
 f=3.263212958871059222

errX = 0.1267

errF = 0.6188

Iteration 9:

 x=2.428849558717576240
 x=-1.655875174933654881
 x=-1.259971976090513479

 f=0.151894835243015364
 f=-0.019188038778540594
 f=0.599601254706248588

errX = 0.0136

errF = 0.0585

Iteration 10:

 x=2.508415279643811502
 x=-1.684939115473996329
 x=-1.133234422365223759

 f=0.019235905449000876
 f=-0.005790633473550177
 f=0.054995347085114506

errX = 1.4954e-04

errF = 6.2642e-04

Iteration 11:

 x=2.514932186684921689
 x=-1.687755323557607712
 x=-1.119595610753047144

 f=0.000202203265315156
 f=-0.000041224669044668
 f=0.000591456437218156

errX = 1.7869e-08

errF = 7.4989e-08

Iteration 12:

 x=2.515002439026511549
 x=-1.687784589308885641
 x=-1.119446066419630181

 f=0.000000024386426293
 f=-0.000000005029495870
 f=0.000000070734383263

errX = 2.2204e-16

errF = 2.2506e-15

Iteration 13:

 x=2.515002447346428749
 x=-1.687784592789972171
 x=-1.119446048550189809

 f=0.000000000000000926
 f=-0.000000000000001658
 f=0.000000000000001208

[amr makhlouf]

the correct answer is obtained using gams software is

28 VARIABLE x.L = 1.253

VARIABLE y.L = 1.166

VARIABLE z.L = 0.278

which is not the same as the result shown above i even have managed to coded it using mathematica but it gives a different results from those one you gents have provided may i know why ?

[Torsten]

The reason is simple: Your function has several roots, like e.g. x^2 - 2 = 0 has roots at + and - sqrt(2).

If you define constraints on the solution (e.g. that all its components should be >= 0), you might be able to make the solution unique.

[amr makhlouf]

I know how to do that using gams can you tell me how to add this constraint using MATLAB to just compare the results ??

[Torsten]

I set xi to the gams solution and your code confirmed it.

Now you can try to set the initial guess in gams to

[2.515002447346428749 ;-1.687784592789972171 ;-1.119446048550189809]

and see whether gams also confirms your MATLAB solution.

syms  x y z  
errX = 10;
errF = 10;
imax = 25;
TolX = 1e-8;
TolF = 1e-8;
f = [ x^2+x-y^2+x*y+z^2-3
    5*(y^2)+y+y*z-x^2+z-7
    x+y+z+y*z-z^3-3];
J = jacobian(f,[x,y,z]);
Jinv = inv(J);
xi = [1.253;1.166;0.278] ;
i = 0;

while (errX > TolX  || errF > TolF) &&   i < imax    % check the condition 
    fi = double(subs(f,[x y z],[xi(1) xi(2) xi(3)]));
    Ji = double(subs(Jinv,[x y z],[xi(1) xi(2) xi(3)]));
    xip1 = xi - Ji*fi;
    i = i + 1;
    errX = max(abs(xip1-xi))
    errF = norm(fi)
    fprintf('Iteration %d:\n',i); 
    fprintf(' x=%.18f\n',xi)
    fprintf(' f=%.18f\n',fi)
    xi = xip1;
end

errX = 3.4777e-04

errF = 0.0044

Iteration 1:

 x=1.252999999999999892
 x=1.165999999999999925
 x=0.278000000000000025

 f=0.001735000000000000
 f=-0.004081000000000000
 f=-0.000336952000000000

errX = 1.7395e-08

errF = 1.5734e-07

Iteration 2:

 x=1.252652227011884722
 x=1.166212281033995657
 x=0.278213752942726544

 f=0.000000047747524948
 f=0.000000149745829800
 f=0.000000007260201776

errX = 4.9960e-16

errF = 8.2827e-16

Iteration 3:

 x=1.252652210852672754
 x=1.166212263638643831
 x=0.278213769040932368

 f=0.000000000000000486
 f=0.000000000000000388
 f=-0.000000000000000547