About solving a system of equations by 'For' loop

8 Apr 2019
1 Answer
Answer Accepted
Updated 13 Apr 2019
4 Views (30 days)

0 votes

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 Comments
Show 2 older comments Hide 2 older comments

Walter Roberson
Walter Roberson on 8 Apr 2019
How many columns do the J variables have?
Mojtaba Mohareri
Mojtaba Mohareri on 8 Apr 2019
J1, J2,J3,J4 are 4*4, 4*1, 1*4* ,1*1 matrices, respectively.
Walter Roberson
Walter Roberson on 8 Apr 2019
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 on 9 Apr 2019
Edited: Mojtaba Mohareri on 9 Apr 2019
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;

Sign in to comment.

Sign in to answer this question.

 Accepted Answer

Walter Roberson
Walter Roberson on 8 Apr 2019

0 votes

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 Comments
Show 11 older comments Hide 11 older comments

Mojtaba Mohareri
Mojtaba Mohareri on 8 Apr 2019
Could you explain to me what to do? I don't really understand it.
Walter Roberson
Walter Roberson on 8 Apr 2019
The first expression of your equations matrix uses J1(1,:) and you say that J1 is a 4x4 array. J1(1,:) would expand to a 1x4 vector. Your first expression is vector valued and so defines 4 equations all of which need to hold at the same time. Is that your intention?
Mojtaba Mohareri
Mojtaba Mohareri on 8 Apr 2019
J1 is a row vector, for example J1=[1,1,1,1] and k1 is also a column vector [k11;k12;k13;k14]. So J1(1,:)*k1 is k11+k12+k13+k14 and is not vector valued.
Mojtaba Mohareri
Mojtaba Mohareri on 8 Apr 2019
I'm so sorry, your are right. J1 a is 4*4 matrix, but J1(1,:) is a vector.
Mojtaba Mohareri
Mojtaba Mohareri on 8 Apr 2019
So, as I told you, J1(1,:)*k1 is a scaler.
Mojtaba Mohareri
Mojtaba Mohareri on 11 Apr 2019
Dear Walter Roberson. Can I solve my problem with other code? I want just to obtain y0 and z0 after 1000 steps and it's not important for me which code I use. Unfortunately I don't know matlab properly. Could you help me to get the answer?
Walter Roberson
Walter Roberson on 11 Apr 2019
It appears that the double() you are using for efficiency is leading to enough loss of precision that you are getting Inf solutions early on, which lead to NaN in further calculation, which lead to your eqns in the C block all immediately evaluating to FALSE before even reaching the solve(), which leads to the solve() returning empty, which leads to subscript out of range.
When I remove those double(), the calculations are much much much much slower, and iterating i to 1000 would probably not be practical.
I will see what I can come up with.
Walter Roberson
Walter Roberson on 11 Apr 2019
When I change the double() to vpa() with 16 digits (equivalent precision but larger range) then by i = 10, the A matrix is composed of values that have magnitude greater than 10^one million and the symbolic numeric solver starts giving up. That is, your equations are overflowing, and it is not due to round-off (same problem happens with 32 digits for example.)
Everything looks fine in the A matrix until i = 6, at which point you start seeing values > 10 after having been no more than +/- 2 before that. i = 7 gives you A matrix values about 1E65. i = 8 gives you A matrix values about 1E8000 . If you were working in double precision, that would be the point at which the values would overflow to infinity and everything after that would start to fail.
You will need to go back and re-check the equations. With the way they are now, you cannot get through 10 iterations without overflowing symbolic numeric range, so there is no way you could get through 1000 iterations.
Your f1 involves z1^4 and the values that result from f1 generally feed back to become the arguments to f1. When z1 becomes greater than 1, that leads to a chain of 4th powers in increases in range, leading to rapid overflow.
Mojtaba Mohareri
Mojtaba Mohareri on 12 Apr 2019
Thak you for your exact explanation and your help. As you told me, I get the answer by i=6 but not for values i greater than 6. Yes, I'll re-check my equations. I think I can remove somehow z1 from the equations but I have to check them again.
Thank you very much for taking the trouble to help me. I really appreciate your help.
Mojtaba Mohareri
Mojtaba Mohareri on 12 Apr 2019
Dear Walter Roberson, I simplify my problem by assuming z0 to be constant 10 during For Loop. I've also changed the double() to vpa() as you told. Now, I get the answer by i=200, but for i=1000. there's an error as
Error using subsref
Index exceeds matrix dimensions. Error in sym/subsref (line 771)
R_tilde = builtin('subsref',L_tilde,Idx);
Error in ROW (line 46)
k2=[B(1);B(2);B(3);B(4)];
Could you tell me what to do for it?
Walter Roberson
Walter Roberson on 12 Apr 2019
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 on 13 Apr 2019
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.

Sign in to comment.

More Answers (0)

Categories

Find more on Numeric Solvers in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!