Using a loop with changing intial value

1 view (last 30 days)
Will Jeter
Will Jeter on 23 Nov 2020
Commented: David Hill on 23 Nov 2020
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

Answers (1)

David Hill
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
  4 Comments
Will Jeter
Will Jeter on 23 Nov 2020
yeah i did and all i get is this
Unable to perform assignment because brace indexing is not
supported for variables of this type.
Error in yrhtd (line 12)
[V{k}, CT{k}] = ode45(@fn, Vspan, IC(k,:));
David Hill
David Hill on 23 Nov 2020
Make sure your workspace is clear. It should work.
clear;

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!