About solving a system of equations by 'For' loop

5 次查看（过去 30 天）
Mojtaba Mohareri
Mojtaba Mohareri 2019-4-8
评论： Mojtaba Mohareri 2019-4-13
I've written 4 equations within a for loop for 10 steps. But there's the following error:
Attempted to access C(1); index out of bounds because numel(C)=0.
Error in ROW (line 59)
k3=[C(1);C(2);C(3);C(4)];
When I use "double()" for any equation, the answer is not satisfactory for me. I was wondering if you could help me about my problem.
In the following, J1,J2,J3,J4, f1,f2,f3,f4 are knowns.
My currenct script looks like this:
h=1/1000;
y0=[2;2;1;0];
z0=10;
for i=1:10
syms k11 k12 k13 k14 l11
k1=[k11;k12;k13;k14];
eqns = [h*(f1(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(1,:)*k1 +J2(1,:)*l11))-k11==0,
h*(f2(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(2,:)*k1 +J2(2,:)*l11))-k12==0,
h*(f3(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(3,:)*k1 +J2(3,:)*l11))-k13 ==0 ,
h*(f4(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(4,:)*k1 +J2(4,:)*l11))-k14 ==0,
g(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J3*k1 +J4*l11)==0];
vars = [k11, k12, k13, k14, l11];
[solk11,solk12,solk13,solk14,soll11] = solve(eqns,vars);
A =double( [solk11,solk12,solk13,solk14,soll11]);
k1=[A(1);A(2);A(3);A(4)];
l11=A(5);
syms k21 k22 k23 k24 l21
k2=[k21;k22;k23;k24];
eqns = [h*(f1(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0+E(5)*l11)+E(11)*(J1(1,:)*k1 +J2(1,:)*l11)+0.44*(J1(1,:)*k2 +J2(1,:)*l21))-k21==0,
h*(f2(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0+E(5)*l11)+E(11)*(J1(2,:)*k1 +J2(2,:)*l11)+0.44*(J1(2,:)*k2 +J2(2,:)*l21))-k22==0,
h*(f3(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0+E(5)*l11)+E(11)*(J1(3,:)*k1 +J2(3,:)*l11)+0.44*(J1(3,:)*k2 +J2(3,:)*l21))-k23 ==0 ,
h*(f4(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0+E(5)*l11)+E(11)*(J1(4,:)*k1 +J2(4,:)*l11)+0.44*(J1(4,:)*k2 +J2(4,:)*l21))-k24 ==0,
g(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0+E(5)*l11)+E(11)*(J3*k1 +J4*l11)+0.44*(J3*k2 +J4*l21)==0];
vars = [k21, k22, k23, k24, l21];
[solk21,solk22,solk23,solk24,soll21] = solve(eqns,vars);
B =double( [solk21,solk22,solk23,solk24,soll21]);
k2=[B(1);B(2);B(3);B(4)];
l21=B(5);
syms k31 k32 k33 k34 l31
k3=[k31;k32;k33;k34];
eqns = [h*(f1(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0+E(6)*l11+E(7)*l21)+E(12)*(J1(1,:)*k1 +J2(1,:)*l11)+E(13)*(J1(1,:)*k2 +J2(1,:)*l21)+0.44*(J1(1,:)*k3 +J2(1,:)*l31))-k31==0,
h*(f2(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0+E(6)*l11+E(7)*l21)+E(12)*(J1(2,:)*k1 +J2(2,:)*l11)+E(13)*(J1(2,:)*k2 +J2(2,:)*l21)+0.44*(J1(2,:)*k3 +J2(2,:)*l31))-k32==0,
h*(f3(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0+E(6)*l11+E(7)*l21)+E(12)*(J1(3,:)*k1 +J2(3,:)*l11)+E(13)*(J1(3,:)*k2 +J2(3,:)*l21)+0.44*(J1(3,:)*k3 +J2(3,:)*l31))-k33==0 ,
h*(f4(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0+E(6)*l11+E(7)*l21)+E(12)*(J1(4,:)*k1 +J2(4,:)*l11)+E(13)*(J1(4,:)*k2 +J2(4,:)*l21)+0.44*(J1(4,:)*k3 +J2(4,:)*l31))-k34 ==0,
g(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0+E(6)*l11+E(7)*l21)+E(11)*(J3*k1 +J4*l11)+E(13)*(J3*k2 +J4*l21)+0.44*(J3*k3 +J4*l31)==0];
vars = [k31, k32, k33, k34, l31];
[solk31,solk32,solk33,solk34,soll31] = solve(eqns,vars);
C =double([solk31,solk32,solk33,solk34,soll31]);
k3=[C(1);C(2);C(3);C(4)];
l31=C(5);
syms k41 k42 k43 k44 l41
k4=[k41;k42;k43;k44];
eqns = [h*(f1(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0+E(8)*l11+E(9)*l21+E(10)*l31)+E(14)*(J1(1,:)*k1 +J2(1,:)*l11)+E(15)*(J1(1,:)*k2 +J2(1,:)*l21)+E(16)*(J1(1,:)*k3 +J2(1,:)*l31)+0.44*(J1(1,:)*k4 +J2(1,:)*l41))-k41==0,
h*(f2(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0+E(8)*l11+E(9)*l21+E(10)*l31)+E(14)*(J1(2,:)*k1 +J2(2,:)*l11)+E(15)*(J1(2,:)*k2 +J2(2,:)*l21)+E(16)*(J1(2,:)*k2 +J2(2,:)*l21)+0.44*(J1(2,:)*k4 +J2(2,:)*l41))-k42==0,
h*(f3(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0+E(8)*l11+E(9)*l21+E(10)*l31)+E(14)*(J1(3,:)*k1 +J2(3,:)*l11)+E(15)*(J1(3,:)*k2 +J2(3,:)*l21)+E(16)*(J1(3,:)*k3 +J2(3,:)*l31)+0.44*(J1(3,:)*k4 +J2(3,:)*l41))-k43 ==0 ,
h*(f4(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0+E(8)*l11+E(9)*l21+E(10)*l31)+E(14)*(J1(4,:)*k1 +J2(4,:)*l11)+E(15)*(J1(4,:)*k2 +J2(4,:)*l21)+E(16)*(J1(4,:)*k3 +J2(4,:)*l31)+0.44*(J1(4,:)*k4 +J2(4,:)*l41))-k44 ==0,
g(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0+E(8)*l11+E(9)*l21+E(10)*l31)+E(14)*(J3*k1 +J4*l11)+E(15)*(J3*k2 +J4*l21)+E(16)*(J3*k3 +J4*l31)+0.44*(J3*k4 +J4*l41)==0];
vars = [k41, k42, k43, k44, l41];
[solk41,solk42,solk43,solk44,soll41] = solve(eqns,vars);
D =double( [solk41,solk42,solk43,solk44,soll41]);
k4=[D(1);D(2);D(3);D(4)];
l41=D(5);
y0=y0+E(1)*k1+E(2)*k2+E(3)*k3+E(4)*k4;
z0=z0+E(1)*l11+E(2)*l21+E(3)*l31+E(4)*l41;
end
  4 个评论
显示 2更早的评论
Walter Roberson
Walter Roberson 2019-4-8
It appears that I will need f1 function, and J* values, and E values, and possibly some others, in order to run the code to test.
Mojtaba Mohareri
Mojtaba Mohareri 2019-4-9
编辑：Mojtaba Mohareri 2019-4-9
Thank you for your help. Yes, of course.
E=[0.3285609536316354, -0.5785609536316354, 0.25, 1, 2.218787467653286, 0, 0, 1.208587690772214, 0.07511610241919324, 0.5, -2.218787467653286, -0.09461966143940745, -0.007913526735718213, -1.870323744195384, -0.09624340112825115 , 0.2726301276675511 ];
f1(y1,y2,y3,y4,z1)=(-(z1^3)/(y3^2))*(3*(y2-y1+y3^(-1)-z1/10)^2+(1/5)*(y2-y1+y3^(-1)-(z1)/10))-y4;
f2(y1,y2,y3,y4,z1)=1/10*z1-y4;
f3(y1,y2,y3,y4,z1)=(z1^3)*(3*(y2-y1+y3^(-1)-z1/10)^2+(1/5)*(y2-y1+y3^(-1)-(z1)/10));
f4(y1,y2,y3,y4,z1)=y1-y3^(-1);
g(y1,y2,y3,y4,z1)=(y1-y3^(-1))^2+y4^2-1/10*z1;
J1(y1,y2,y3,y4,z1)=jacobian([f1,f2,f3,f4],[y1,y2,y3,y4]);
J1=J1(2,2,1,0,10);
J2(y1,y2,y3,y4,z1)=jacobian([f1,f2,f3,f4],[z1]);
J2=J2(2,2,1,0,10);
J3(y1,y2,y3,y4,z1)=jacobian([g],[y1,y2,y3,y4]);
J3=J3(2,2,1,0,10);
J4(y1,y2,y3,y4,z1)=jacobian([g],[z1]);
J4=J4(2,2,1,0,10);
y0=[2;2;1;0];
z0=10;

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2019-4-8
J1(1,:) is a vector of 4 elements. J2(1,:) is a scalar. When you use both in a single == expression then you are defining 4 simultaneous equations that vary only in J1 values, but you expect that single values of each of the variables will be able to satisfy all of those 4 equations simultaneously. That can only work if the defined equations are degenerate, where the J1 term is being multiplied by 0.
  13 个评论
显示 11更早的评论
Walter Roberson
Walter Roberson 2019-4-12
Please post your current code. When I try your code as posted before except with the change to z0 commented out, then it stops at i = 102 with "integer too large in context" because it has reached values that exceed 10 to the three and a half million power.
Mojtaba Mohareri
Mojtaba Mohareri 2019-4-13
Yes, you are right, but I don't know what the problem is. In your opinion, Is there something wrong with my equations? or with my codes? I re-check my equations every time:(
Yes, sure. My current code:
syms b1 b2 b_3 b_4 a21 a31 a32 a41 a42 a43 g21 g31 g32 g41 g42 g43 w12 w22 y1 y2 y3 y4 z1 k11 k12 k13 k14 k21 k22 k23 k24 k31 k32 k33 k34 k41 k42 k43 k44 l11 l21 l31 l41 x
k=[k11;k12;k13;k14];
a=0;
h=1/1000;
E=[0.3285609536316354, -0.5785609536316354, 0.25, 1, 2.218787467653286, 0, 0, 1.208587690772214, 0.07511610241919324, 0.5, -2.218787467653286, -0.09461966143940745, -0.007913526735718213, -1.870323744195384, -0.09624340112825115 , 0.2726301276675511 ];
f1(y1,y2,y3,y4,z1)=(-(z1^3)/(y3^2))*(3*(y2-y1+y3^(-1)-z1/10)^2+(1/5)*(y2-y1+y3^(-1)-(z1)/10))-y4;
f2(y1,y2,y3,y4,z1)=1/10*z1-y4;
f3(y1,y2,y3,y4,z1)=(z1^3)*(3*(y2-y1+y3^(-1)-(z1)/10)^2+(1/5)*(y2-y1+y3^(-1)-(z1)/10));
f4(y1,y2,y3,y4,z1)=y1-y3^(-1);
g(y1,y2,y3,y4,z1)=(y1-y3^(-1))^2+y4^2-1/10*z1;
J1(y1,y2,y3,y4,z1)=jacobian([f1,f2,f3,f4],[y1,y2,y3,y4]);
J1=J1(2,2,1,0,10);
J2(y1,y2,y3,y4,z1)=jacobian([f1,f2,f3,f4],[z1]);
J2=J2(2,2,1,0,10);
J3(y1,y2,y3,y4,z1)=jacobian([g],[y1,y2,y3,y4]);
J3=J3(2,2,1,0,10);
J4=-1/10;
y0=[2;2;1;0];
z0=10;
for i=1:500
syms k11 k12 k13 k14 l11
k1=[k11;k12;k13;k14];
eqns = [h*(f1(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(1,:)*k1 +J2(1,:)*l11))-k11==0,
h*(f2(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(2,:)*k1 +J2(2,:)*l11))-k12==0,
h*(f3(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(3,:)*k1 +J2(3,:)*l11))-k13 ==0 ,
h*(f4(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J1(4,:)*k1 +J2(4,:)*l11))-k14 ==0,
g(y0(1),y0(2),y0(3),y0(4),z0)+0.44*(J3*k1 +J4*l11)==0];
vars = [k11, k12, k13, k14, l11];
[solk11,solk12,solk13,solk14,soll11] = solve(eqns,vars);
A =vpa([solk11,solk12,solk13,solk14,soll11]);
k1=[A(1);A(2);A(3);A(4)];
l11=A(5);
syms k21 k22 k23 k24 l21
k2=[k21;k22;k23;k24];
eqns = [h*(f1(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0)+E(11)*(J1(1,:)*k1 +J2(1,:)*l11)+0.44*(J1(1,:)*k2 +J2(1,:)*l21))-k21==0,
h*(f2(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0)+E(11)*(J1(2,:)*k1 +J2(2,:)*l11)+0.44*(J1(2,:)*k2 +J2(2,:)*l21))-k22==0,
h*(f3(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0)+E(11)*(J1(3,:)*k1 +J2(3,:)*l11)+0.44*(J1(3,:)*k2 +J2(3,:)*l21))-k23 ==0 ,
h*(f4(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0)+E(11)*(J1(4,:)*k1 +J2(4,:)*l11)+0.44*(J1(4,:)*k2 +J2(4,:)*l21))-k24 ==0,
g(y0(1)+E(5)*k1(1),y0(2)+E(5)*k1(2),y0(3)+E(5)*k1(3),y0(4)+E(5)*k1(4),z0)+E(11)*(J3*k1 +J4*l11)+0.44*(J3*k2 +J4*l21)==0];
vars = [k21, k22, k23, k24, l21];
[solk21,solk22,solk23,solk24,soll21] = solve(eqns,vars);
B =vpa( [solk21,solk22,solk23,solk24,soll21]);
k2=[B(1);B(2);B(3);B(4)];
l21=B(5);
syms k31 k32 k33 k34 l31
k3=[k31;k32;k33;k34];
eqns = [h*(f1(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0)+E(12)*(J1(1,:)*k1 +J2(1,:)*l11)+E(13)*(J1(1,:)*k2 +J2(1,:)*l21)+0.44*(J1(1,:)*k3 +J2(1,:)*l31))-k31==0,
h*(f2(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0)+E(12)*(J1(2,:)*k1 +J2(2,:)*l11)+E(13)*(J1(2,:)*k2 +J2(2,:)*l21)+0.44*(J1(2,:)*k3 +J2(2,:)*l31))-k32==0,
h*(f3(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0)+E(12)*(J1(3,:)*k1 +J2(3,:)*l11)+E(13)*(J1(3,:)*k2 +J2(3,:)*l21)+0.44*(J1(3,:)*k3 +J2(3,:)*l31))-k33 ==0 ,
h*(f4(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0)+E(12)*(J1(4,:)*k1 +J2(4,:)*l11)+E(13)*(J1(4,:)*k2 +J2(4,:)*l21)+0.44*(J1(4,:)*k3 +J2(4,:)*l31))-k34 ==0,
g(y0(1)+E(6)*k1(1)+E(7)*k2(1),y0(2)+E(6)*k1(2)+E(7)*k2(2),y0(3)+E(6)*k1(3)+E(7)*k2(3),y0(4)+E(6)*k1(4)+E(7)*k2(4),z0)+E(11)*(J3*k1 +J4*l11)+E(13)*(J3*k2 +J4*l21)+0.44*(J3*k3 +J4*l31)==0];
vars = [k31, k32, k33, k34, l31];
[solk31,solk32,solk33,solk34,soll31] = solve(eqns,vars);
C =vpa([solk31,solk32,solk33,solk34,soll31]);
k3=[C(1);C(2);C(3);C(4)];
l31=C(5);
syms k41 k42 k43 k44 l41
k4=[k41;k42;k43;k44];
eqns = [h*(f1(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0)+E(14)*(J1(1,:)*k1 +J2(1,:)*l11)+E(15)*(J1(1,:)*k2 +J2(1,:)*l21)+E(16)*(J1(1,:)*k3 +J2(1,:)*l31)+0.44*(J1(1,:)*k4 +J2(1,:)*l41))-k41==0,
h*(f2(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0)+E(14)*(J1(2,:)*k1 +J2(2,:)*l11)+E(15)*(J1(2,:)*k2 +J2(2,:)*l21)+E(16)*(J1(2,:)*k2 +J2(2,:)*l21)+0.44*(J1(2,:)*k4 +J2(2,:)*l41))-k42==0,
h*(f3(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0)+E(14)*(J1(3,:)*k1 +J2(3,:)*l11)+E(15)*(J1(3,:)*k2 +J2(3,:)*l21)+E(16)*(J1(3,:)*k3 +J2(3,:)*l31)+0.44*(J1(3,:)*k4 +J2(3,:)*l41))-k43 ==0 ,
h*(f4(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0)+E(14)*(J1(4,:)*k1 +J2(4,:)*l11)+E(15)*(J1(4,:)*k2 +J2(4,:)*l21)+E(16)*(J1(4,:)*k3 +J2(4,:)*l31)+0.44*(J1(4,:)*k4 +J2(4,:)*l41))-k44 ==0,
g(y0(1)+E(8)*k1(1)+E(8)*k2(1)+E(10)*k3(1),y0(2)+E(8)*k1(2)+E(8)*k2(2)+E(10)*k3(2),y0(3)+E(8)*k1(3)+E(8)*k2(3)+E(10)*k3(3),y0(4)+E(8)*k1(4)+E(8)*k2(4)+E(10)*k3(4),z0)+E(14)*(J3*k1 +J4*l11)+E(15)*(J3*k2 +J4*l21)+E(16)*(J3*k3 +J4*l31)+0.44*(J3*k4 +J4*l41)==0];
vars = [k41, k42, k43, k44, l41];
[solk41,solk42,solk43,solk44,soll41] = solve(eqns,vars);
D =vpa( [solk41,solk42,solk43,solk44,soll41]);
k4=[D(1);D(2);D(3);D(4)];
l41=D(5);
y0=y0+E(1)*k1+E(2)*k2+E(3)*k3+E(4)*k4;
y{i}=y0;
w(i)=y0(4);
z0=10;
end
Thank you very much.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Numeric Solvers 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by