Cant get my loop to work
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Will Jeter
am 23 Nov. 2020
Kommentiert: Jan
am 24 Nov. 2020
apparently there's something wrong in line 12 and i keep getting a debugging page for ode45 when i try to run it
% Initial conditions
Ca0 = 2;
Cb0 = 4;
Cc0 = 0;
Cd0 = 0;
T0 = 300:100:800;
IC = [repmat([Ca0 Cb0 Cc0 Cd0],length(T0),1),T0'];
% V range
Vspan = [0 10];
% Call ode solver
for k=1:size(IC,1)
[V{k}, CT{k}] = ode45(@fn, Vspan, IC(k,:));
figure;
plot(V{k},CT{k}(:,1),V{k},CT{k}(:,2),V{k},CT{k}(:,3),V{k},CT{k}(:,4));
xlabel('V (dm^3)'),ylabel('Ca.Cb.Cc,Cd (mol/dm^3)');
legend('Ca','Cb','Cc','Cd');
figure
plot(V{k},CT{k}(:,5));
grid on;
xlabel('V(dm^3)'),ylabel('T(K)');
end
% ODE function
function dCTdV = fn(~,CT)
% Data
Cpa = 20;
dh1a = 20000;
dh2a = -10000;
vo = 10;
% Extract parameters
Ca = CT(1);
Cb = CT(2);
Cc = CT(3);
Cd = CT(4);
T = CT(5);
% k and r functions
k1a = 0.001*exp(-5000/1.987*(1/T - 1/300));
k2a = 0.001*exp(-7500/1.987*(1/T - 1/300));
r1a = -k1a*Ca*Cb^2;
r2a = -k2a*Ca*Cc;
% odes
dCTdV = [(r1a+r2a)/vo;
2*r1a/vo;
(-2*r1a+r2a)/vo;
-2*r2a/vo;
(r1a*dh1a + r2a*dh2a)/((Ca+Cb+3*Cc+4*Cd)*vo*Cpa)];
end
2 Kommentare
Jan
am 24 Nov. 2020
I recommend to mention, why you assume, that there is a problem in line 12. Do you get an error message? Where does the debugger stop?
Akzeptierte Antwort
Alan Stevens
am 24 Nov. 2020
Replace your k loop by
for k=1:size(IC,1)
[V, CT] = ode45(@fn, Vspan, IC(k,:));
figure;
plot(V,CT(:,1),V,CT(:,2),V,CT(:,3),V,CT(:,4));
xlabel('V (dm^3)'),ylabel('Ca.Cb.Cc,Cd (mol/dm^3)');
legend('Ca','Cb','Cc','Cd');
figure
plot(V,CT(:,5));
grid on;
xlabel('V(dm^3)'),ylabel('T(K)');
end
0 Kommentare
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!