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

1 回表示 (過去 30 日間)
Akshay Pratap Singh
Akshay Pratap Singh 2019 年 1 月 23 日
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)])

採用された回答

Stephan
Stephan 2019 年 1 月 23 日
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

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeLinear Algebra についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by