Using a loop with changing intial value

조회 수: 6 (최근 30일)
Will Jeter
Will Jeter 2020년 11월 23일
댓글: David Hill 2020년 11월 23일
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

답변 (1개)

David Hill
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
  댓글 수: 4
Will Jeter
Will Jeter 2020년 11월 23일
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 2020년 11월 23일
Make sure your workspace is clear. It should work.
clear;

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Creating and Concatenating Matrices에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by