Trying to code from Polymath

Question

Will Jeter 2020-11-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/658048-trying-to-code-from-polymath

回答： Alan Stevens 2020-11-23

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Stevens 2020-11-23

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/658048-trying-to-code-from-polymath#answer_553098

在 MATLAB Online 中打开

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