How to find the jacobian of the system of nonlinear equations?
2 views (last 30 days)
Show older comments
Akshay Pratap Singh
on 23 Jan 2019
Commented: Akshay Pratap Singh
on 23 Jan 2019
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)])
0 Comments
Accepted Answer
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)
See Also
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!