Using a loop with changing intial value
1 回表示 (過去 30 日間)
古いコメントを表示
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 件のコメント
回答 (1 件)
David Hill
2020 年 11 月 23 日
% 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
参考
カテゴリ
Help Center および File Exchange で Ordinary Differential Equations についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!