Using a loop with changing intial value
1 view (last 30 days)
Show older comments
I need to use a loop to be used to solve the ODE with various T0 which needs to range from 300 to 800. Right now i only know how to do one T0 at a time, but i need to evalute all from 300 to 800. Please help.
% Initial conditions
Ca0 = 2;
Cb0 = 4;
Cc0 = 0;
Cd0 = 0;
T0 = 600;
IC = [Ca0 Cb0 Cc0 Cd0 T0];
% V range
Vspan = [0 10];
% Call ode solver
[V, CT] = ode45(@fn, Vspan, IC);
% Extract parameters
Ca = CT(:,1);
Cb = CT(:,2);
Cc = CT(:,3);
Cd = CT(:,4);
T = CT(:,5);
% Plot graphs
figure
plot(V,Ca,V,Cb,V,Cc,V,Cd),grid
xlabel('V (dm^3)'),ylabel('Ca.Cb.Cc,Cd (mol/dm^3)')
legend('Ca','Cb','Cc','Cd')
figure
plot(V,T),grid
xlabel('V(dm^3)'),ylabel('T(K)')
% 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
0 Comments
Answers (1)
David Hill
on 23 Nov 2020
% 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
See Also
Categories
Find more on Ordinary Differential Equations 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!