MATLAB: How to find the jacobian of the system of nonlinear equations

Question

I wrote a following program, but it is giving some errors?

clear all
syms x
h=4;
gma=18.4;
ka1=0.2;
kp1=8.76
ka2=0.2;
kp2=8.76;
sindel=0.437;
cosdel=0.9;
pa1=ka1*gma*(h*x(1)+0.5*x(1)^2);
pp1=kp1*gma*0.5*x(1)^2;
pa2=ka2*gma*(x(1)*x(2)+0.5*x(2)^2);
pp2=kp2*(x(2)*(h+x(1))+0.5*x(2)^2);
za1=(0.5*h*x(1)^2+(x(1)^3)/6)/(h*x(1)+0.5*x(1)^2);
zp2=(0.5*(h+x(1))*x(2)^2+((x(2)^3)/3))/((h+x(1))*x(2)+0.5*x(2)^2);
za2=(0.5*x(1)*x(2)^2+((x(2)^3)/3))+(x(1)*x(2)+0.5*x(2)^2);
e1=pp1*sindel-pa1*sindel-pp2*sindel-pa2*sindel;
e2=pp1*cosdel+pa2*cosdel-pa1*cosdel-pp2*cosdel-x(3);
e3=pp1*cosdel*(x(1)/3)+pp2*cosdel*zp2-pa1*cosdel*za1-pa2*cosdel*za2-x(3)*(x(1)+(h/3));
g=[e1; e2; e3];
J=jacobian([e1, e2, e3], [x(1), x(2), x(3)])

Answer 1

Hi,.

try:

syms x1 x2 x3
h=4;
gma=18.4;
ka1=0.2;
kp1=8.76
ka2=0.2;
kp2=8.76;
sindel=0.437;
cosdel=0.9;
pa1=ka1*gma*(h*x1+0.5*x1^2);
pp1=kp1*gma*0.5*x1^2;
pa2=ka2*gma*(x1*x2*+0.5*x2^2);
pp2=kp2*(x2*(h+x1)+0.5*x2^2);
za1=(0.5*h*x1^2+(x1^3)/6)/(h*x1+0.5*x1^2);
zp2=(0.5*(h+x1)*x2^2+((x2^3)/3))/((h+x1)*x2+0.5*x2^2);
za2=(0.5*x1*x2^2+((x2^3)/3))+(x1*x2+0.5*x2^2);
e1=pp1*sindel-pa1*sindel-pp2*sindel-pa2*sindel;
e2=pp1*cosdel+pa2*cosdel-pa1*cosdel-pp2*cosdel-x3;
e3=pp1*cosdel*(x1/3)+pp2*cosdel*zp2-pa1*cosdel*za1-pa2*cosdel*za2-x3*(x1+(h/3));
g=[e1; e2; e3];
J=jacobian([e1, e2, e3], [x1, x2, x3])

results in:

J =
[                                                                                                                                                                                                                                                                                                                                                                                                                                                                     (1075457*x1)/15625 - (95703*x2)/25000 - (10051*x2^3)/12500 - 20102/3125,                                                                                                                                                                                                                                                                                                                                              - (95703*x1)/25000 - (95703*x2)/25000 - (30153*x1*x2^2)/12500 - 95703/6250,          0]
[                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  (207*x2^3)/125 - (1971*x2)/250 + (88596*x1)/625 - 1656/125,                                                                                                                                                                                                                                                                                                                                                            (621*x1*x2^2)/125 - (1971*x2)/250 - (1971*x1)/250 - 3942/125,         -1]
[ (44298*x1^2)/625 - x3 - (207*x2^3*(x1*x2 + (x1*x2^2)/2 + x2^2/2 + x2^3/3))/125 - (1656*x1)/125 + (x2^2*((1971*x2*(x1 + 4))/250 + (1971*x2^2)/500))/(2*(x2^2/2 + (x1 + 4)*x2)) - (((414*x1)/125 + 1656/125)*(x1^3/6 + 2*x1^2))/(x1^2/2 + 4*x1) + (1971*x2*(x2^3/3 + (x1/2 + 2)*x2^2))/(250*(x2^2/2 + (x1 + 4)*x2)) - (207*x1*x2^3*(x2^2/2 + x2))/125 + (((207*x1^2)/125 + (1656*x1)/125)*(x1^3/6 + 2*x1^2)*(x1 + 4))/(x1^2/2 + 4*x1)^2 - (x2*(x2^3/3 + (x1/2 + 2)*x2^2)*((1971*x2*(x1 + 4))/250 + (1971*x2^2)/500))/(x2^2/2 + (x1 + 4)*x2)^2, ((x2^3/3 + (x1/2 + 2)*x2^2)*((1971*x1)/250 + (1971*x2)/250 + 3942/125))/(x2^2/2 + (x1 + 4)*x2) - (207*x1*x2^3*(x1 + x2 + x1*x2 + x2^2))/125 - (621*x1*x2^2*(x1*x2 + (x1*x2^2)/2 + x2^2/2 + x2^3/3))/125 + ((2*x2*(x1/2 + 2) + x2^2)*((1971*x2*(x1 + 4))/250 + (1971*x2^2)/500))/(x2^2/2 + (x1 + 4)*x2) - ((x2^3/3 + (x1/2 + 2)*x2^2)*((1971*x2*(x1 + 4))/250 + (1971*x2^2)/500)*(x1 + x2 + 4))/(x2^2/2 + (x1 + 4)*x2)^2, - x1 - 4/3]

Best regards

Stephan