for loop, while loop

1 view (last 30 days)
Delia Bosshart on 24 May 2021
Commented: Delia Bosshart on 27 May 2021
Is there a way to automatize the following code with a for or while loop? such that I wouldn't have to rewrite the code for each iterating time step? I would like a loop in which it would give me the result (x1 x2 x3 x4 ...) for every time step (t2 t3 t4 t5 ...). I am confused because variable d is dependent on the result of each step (see code below). Also how can I save the data of the output if the varible x is a sym? It only shows me the value for x in the command window but not as a value in the workspace. I would like to have a tabula of all x values as output that I can save.
syms x1 x2 x3
%t2
a = trapz(t(1:2),q_ext(1:2))/t(2); %q_ext,kumuliert
b = trapz(t(1:2),q_losses(1:2))/t(2); %q_losses,kumuliert
c = 0; %betta_char, konstant in m/s
d = (x1-c).*(t(2)-t(1)); %h_char
e = 232.87.*d.^(-0.46); %rho_char
eqn = (a - b + 6.96.*x1 + (6-(31*e)/1000).*x1/60*1000 + 31.*e/1000.*c/60*1000)*0.0086 - 0.1475 - x1 == 0;
sol_x1 = vpasolve(eqn, x1)
%t3
a = trapz(t(1:3),q_ext(1:3))/t(3);
b = trapz(t(1:3),q_losses(1:3))/t(3);
c = 0;
d = (x1-c).*(t(2)-t(1))+(x2-c).*(t(3)-t(2));
e = 232.87.*d.^(-0.46);
eqn = (a - b + 6.96.*x2 + (6-(31*e)/1000).*x2/60*1000 + 31.*e/1000.*c/60*1000)*0.0086 - 0.1475 - x2 == 0;
sol_x2 = vpasolve(eqn, x2)
%t4
a = trapz(t(1:4),q_ext(1:4))/t(4);
b = trapz(t(1:4),q_losses(1:4))/t(4);
c = 0;
d = (x1-c).*(t(2)-t(1))+(x2-c).*(t(3)-t(2))+(x3-c).*(t(4)-t(3));
e = 232.87.*d.^(-0.46);
eqn = (a - b + 6.96.*x3 + (6-(31*e)/1000).*x3/60*1000 + 31.*e/1000.*c/60*1000)*0.0086 - 0.1475 - x3 == 0;
sol_x3 = vpasolve(eqn, x3)

David Hill on 24 May 2021
Something like this should work. Could not check, since all variables were not provided.
syms x1 x2 x3
x=[x1 x2 x3];
for k=1:3
a = trapz(t(1:k+1),q_ext(1:k+1))/t(k+1);
b = trapz(t(1:k+1),q_losses(1:k+1))/t(k+1);
c = 0;
d = sum(diff(t(1:k+1)).*(x(1:k)-c));
e = 232.87*d^(-0.46);
eqn = (a - b + 6.96*x(k) + (6-(31*e)/1000)*x(k)/60*1000 + 31*e/1000*c/60*1000)*0.0086 - 0.1475 - x(k) == 0;
sol{k} = vpasolve(eqn, x(k));
end
Delia Bosshart on 27 May 2021
Hello David. When using the code, it doesn't show the values for x but only indicates them as syms:
sol =
1×3 cell array
{1×1 sym} {1×1 sym} {1×1 sym}
Why is that? and how can I display the wanted values?
t = [0;12;18;25];
q_ext = [30;50;83;132];
syms x1 x2 x3
x=[x1 x2 x3];
for k=1:3
a = trapz(t(k:k+1),q_ext(k:k+1))/(t(k+1)-t(k));
b = 0;
c = 0;
d = sum(diff(t(k:k+1)).*(x(k)));
e = 232.87*d^(-0.46);
eqn = (a + b + 6.96*x(k) + (6-(31*e)/1000)*x(k)/60*1000 + 31*e/1000*c/60*1000)*0.0086 - 0.1475 - x(k) == 0;
sol{k} = vpasolve(eqn, x(k))
end