How to find the jacobian of the system of nonlinear equations?

1 view (last 30 days)
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)])

Accepted Answer

Stephan
Stephan on 23 Jan 2019
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

More Answers (0)

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!