solve multipe eq in a loop
Show older comments
hello
would you plz fix this code
i dont know how to solve multiple variables in a loop
clc
clear all
dt=1 ;V=1000*1000*100 ; pro=0.15 ; Co=10^-6 ; alpha =5.615; Bo=1.35;
beta =1.127 ;Ax=1000*100 ; mio=1 ;Dx=1000;Kx=20*10^-3;
cte=(V*pro*Co)/(alpha*Bo*dt);
T=(beta * Ax* Kx)/(mio *Bo*Dx);
p1(1:360)=5000 ; p2(1)=5000; p3(1)=5000;p4(1)=5000 ; p5(1)=5000; p6(1)=5000;
p7(1)=5000 ; p8(1)=5000; p9(1)=5000;
syms p2(i+1) p3(i+1) p4(i+1) p5(i+1) p6(i+1) p7(i+1) p8(i+1) p9(i+1)
for i=1:1:359
eqn = [
p3(i+1)*T-(cte +2*T)*p2(i+1)+T*p1(i+1)==-(cte *p2(i));
p4(i+1)*T-(cte +2*T)*p3(i+1)+T*p2(i+1)==-(cte *p2(i)-250);
p5(i+1)*T-(cte +2*T)*p4(i+1)+T*p3(i+1)==-(cte *p3(i));
p6(i+1)*T-(cte +2*T)*p5(i+1)+T*p4(i+1)==-(cte *p4(i));
p7(i+1)*T-(cte +2*T)*p6(i+1)+T*p5(i+1)==-(cte *p5(i));
p8(i+1)*T-(cte +2*T)*p7(i+1)+T*p6(i+1)==-(cte *p6(i)-250);
p9(i+1)*T-(cte +2*T)*p8(i+1)+T*p7(i+1)==-(cte *p7(i));
-(cte +2*T)*p9(i+1)+T*p8(i+1)==-(cte *p8(i))
];
S = solve(eqn, [p2(i+1), p3(i+1), p4(i+1), p5(i+1), p6(i+1), p7(i+1),
p8(i+1),p9(i+1)]);
p2(i+1) = S.p2(i+1);
p3(i+1) = S.p3(i+1);
p4(i+1) = S.p4(i+1);
p5(i+1) = S.p5(i+1);
p6(i+1) = S.p6(i+1);
p7(i+1) = S.p7(i+1);
p8(i+1) = S.p8(i+1);
p9(i+1) = S.p9(i+1);
end
time_d = {'90';'180';'270';'360'};
p1= [p1(90);p1(180);p1(270);p1(360)];
p2 = [p2(90);p2(180);p2(270);p2(360)];
p3 = [p3(90);p3(180);p3(270);p3(360)];
p4 = [p4(90);p4(180);p4(270);p4(360)];
p5 = [p5(90);p5(180);p5(270);p5(360)];
p6= [p6(90);p6(180);p6(270);p6(360)];
p7= [p7(90);p7(180);p7(270);p7(360)];
p8= [p8(90);p8(180);p8(270);p8(360)];
p9= [p9(90);p9(180);p9(270);p9(360)];
T1 = table(time_d,p1,p2,p3,p4,p5,p6,p7,p8,p9)
1 Comment
KALYAN ACHARJYA
on 1 Jan 2021
Please read the following thread
https://in.mathworks.com/matlabcentral/answers/304528-tutorial-why-variables-should-not-be-named-dynamically-eval
Answers (1)
Alan Stevens
on 1 Jan 2021
Your equations are linear in p (p2 to p9, that is, since p1 is constant for all times), so the following is a simple way to solve for them:
dt=1 ;V=1000*1000*100 ; pro=0.15 ; Co=10^-6 ; alpha =5.615; Bo=1.35;
beta =1.127 ;Ax=1000*100 ; mio=1 ;Dx=1000;Kx=20*10^-3;
cte=(V*pro*Co)/(alpha*Bo*dt);
T=(beta * Ax* Kx)/(mio *Bo*Dx);
N = 360;
p = zeros(8,N);
p1=5000 ; p(:,1) = 5000;
k = -(cte + 2*T);
% Equations are linear in p, so solve using M*p(i) = -cte*p(i-1) - K
% or p(i) = M\(-cte*p(i-1)-K), where M and K are given below.
M = [k, T, 0, 0, 0, 0, 0, 0;
T, k, T, 0, 0, 0, 0, 0;
0, T, k, T, 0, 0, 0, 0;
0, 0, T, k, T, 0, 0, 0;
0, 0, 0, T, k, T, 0, 0;
0, 0, 0, 0, T, k, T, 0;
0, 0, 0, 0, 0, T, k, T;
0, 0, 0, 0, 0, 0, T, k];
K = [T*p1; 250; 0; 0; 0; 250; 0; 0];
for i = 2:N
p(:,i) = M\(-cte*p(:,i-1) - K);
end
t = 1:N;
plot(t,p),grid
xlabel('t'),ylabel('p2 to p9')
legend('p2','p3','p4','p5','p6','p7','p8','p9')
Note that the eight rows of p represent p2 to p9.
Categories
Find more on Code Performance 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!