Why is my code not working how I want it to?

2 次查看（过去 30 天）
Marlon Brutus
Marlon Brutus 2017-11-27
评论： MHZ 2017-11-28
Good News Everyone!
My matlab code is malfunctioning and I have no idea why... So I am trying to convert this Polymath code into a Matlab code and it is not working out. I attach the Polymath code and what is supposed to come out but I am just not getting it. Any help will be greatly appreciated. :D
This is my code:
function nonisothermal
vspan=[0 1]; %volume of the reactor dm3
x0=[0 1035]; %Conversion and Temperature(in K)
[v x]=ode45(@thermalfun,vspan,x0);
subplot(2,1,1)
plot(v,x(:,1),'m--','linewidth',1)
xlabel('Volume of the Reactor')
title('Conversion')
grid on
subplot(2,1,2)
plot(v,x(:,2),'k--','linewidth',1)
xlabel('Volume of the Reactor')
title('Temperature')
grid on
end
function dXdV=thermalfun(v,x)
Cao=18.8;
Fao=0.0376;
Cpa=163;
delCp=-9;
To=1035;
delH=80770+delCp*(x(2)-298);
ra=-Cao*3.58*exp(34222*(1/To-1/x(2)))*(1-x(1))*(To/x(2))/(1+x(1));
Ua=165000;
Ta=1150;
dXdV=[-ra/Fao; (Ua*(Ta-x(2))+ra*delH)/(Fao*(Cpa+x(1)*delCp))]
end
  1 个评论
Star Strider
Star Strider 2017-11-28
I do not see any significant problems with your code, with respect to the published equations. There is an insignificant problem in that ‘Ua’ should be 16500, and that the second plot is ‘-ra’ as a function of ‘v’, so you need to calculate ‘ra’ outside the ode45 call with ‘x(:,1)’, ‘x(:,2)’ and the relevant constants.
Check to see the ODE solver the authors used. The MATLAB numerical solvers are the best I’ve encountered, however some papers use much less sophisticated solvers, leading to results that are difficult to reproduce with the MATLAB solvers. (I invariably trust the MATLAB results.) Also experiment with a much shorter upper limit for ‘vspan’.

请先登录，再进行评论。

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

回答（1 个）

MHZ
MHZ 2017-11-27
编辑：per isakson 2017-11-28
I tried out your code and ran it fine. First make sure to save 2 files in a location in matlab path. Just creat a folder to save in and then go to Matlab home> set path> add folder then hit save and close. Although you do not need to save the first one as a function and can do a normal m-script, save the following as nonisothermal
function nonisothermal
vspan=[0 1]; %volume of the reactor dm3
x0=[0 1035]; %Conversion and Temperature(in K)
[v x]=ode45(@thermalfun,vspan,x0);
subplot(2,1,1)
plot(v,x(:,1),'m--','linewidth',1)
xlabel('Volume of the Reactor')
title('Conversion')
grid on
subplot(2,1,2)
plot(v,x(:,2),'k--','linewidth',1)
xlabel('Volume of the Reactor')
title('Temperature')
grid on
end
Now save the following as thermalfun
function dXdV=thermalfun(v,x)
Cao=18.8;
Fao=0.0376;
Cpa=163;
delCp=-9;
To=1035;
delH=80770+delCp*(x(2)-298);
ra=-Cao*3.58*exp(34222*(1/To-1/x(2)))*(1-x(1))*(To/x(2))/(1+x(1));
Ua=165000;
Ta=1150;
dXdV=[-ra/Fao; (Ua*(Ta-x(2))+ra*delH)/(Fao*(Cpa+x(1)*delCp))]
end
When done, go to the command window and type nonisothermal and hit enter.
  2 个评论
Marlon Brutus
Marlon Brutus 2017-11-28
I don't think you understand... The code is supposed to come with results that are pretty similar to what is attached above.
MHZ
MHZ 2017-11-28
Your problem is in the implementation of the integration equation. You are using x(2) instead of T. According to Matlab, x(2) will be 1035 which is your final and not the value of T. You need to separate dT/dV and actually use T.

请先登录，再进行评论。

类别

在 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!

Translated by