Trying to code from Polymath

24 views (last 30 days)
Will Jeter
Will Jeter on 23 Nov 2020
Answered: Alan Stevens on 23 Nov 2020
Im supposed to turn this Polymath code into a Matlab code but Im not sure how. Please help. Its supposed to use an ode solver and a graph

Accepted Answer

Alan Stevens
Alan Stevens on 23 Nov 2020
Here's a basic MATLAB version.
% Initial conditions
Ca0 = 2;
Cb0 = 4;
Cc0 = 0;
Cd0 = 0;
T0 = 800;
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'),ylabel('Ca.Cb.Cc,Cd')
legend('Ca','Cb','Cc','Cd')
figure
plot(V,T),grid
xlabel('V'),ylabel('T')
% 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
I'll leave you to construct the tabular form of the results.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!