About solving a system of equations by 'For' loop
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Mojtaba Mohareri
am 8 Apr. 2019
Kommentiert: Mojtaba Mohareri
am 13 Apr. 2019
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 Kommentare
Walter Roberson
am 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.
Akzeptierte Antwort
Walter Roberson
am 8 Apr. 2019
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 Kommentare
Walter Roberson
am 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.
Weitere Antworten (0)
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!