Tek-Tips is the largest IT community on the Internet today!

Members share and learn making Tek-Tips Forums the best source of peer-reviewed technical information on the Internet!

  • Congratulations JAE on being selected by the Eng-Tips community for having the most helpful posts in the forums last week. Way to Go!

ODE45 function and graphs not compiling.

Status
Not open for further replies.

anthony0028

Chemical
Joined
Mar 15, 2015
Messages
2
Location
US
I am having issues getting this code to run. I actually received this code back in college when taking a reactor design class, now it will not run. I am receiving several errors which I have failed to resolve. what is wrong with the coding?

I plan to use this as a template for a reactor at the chemical plant I work at. I will ultimately modify the equations, and add an additional differential equation, that is, dx/dw=f(x,p,t), dp/dw=f(x,p,t) , dt/dw=f(x,p,t).

here are the errors when I run the ODE45 script:
>> exampleGas_graphs
[highlight #EF2929]Error using odearguments (line 90)
EXAMPLEGASPD_X_2013 must return a column vector.

Error in ode45 (line 113)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in exampleGas_graphs (line 1)
[w,y]=ode45('ExampleGasPD_X_2013',[0 70],[1 0]);



The function file is as follows:
--------------------------------------
function dy=ExampleGasPD_X_2013(w,y);
dy=zeros(1,0);
alpha=0.01;
vo=10;
Fao=2; Fbo=5; Fto=Fao+Fbo;
T1=(73000*y(2)+53550);
T2=105+30*y(2);
T=(73000*y(2)+53550)/(105+(30*y(2)));
x=y(2);
Cto=Fto/vo;
ko=0.001; Ea=10000; R1=1.987;
To=300;
kprime=ko*exp((Ea/R1)*((1/To)-(1/T)));
Ca=Cto*(1-y(2))*y(1)/(3.5-y(2))*To/T;
Cb=Cto*(2.5-2*y(2))*y(1)/(3.5-y(2))*To/T;
Cc=Cto*(2*y(2))*y(1)/(3.5-y(2))*To/T;
rate=kprime*(Ca*Cb*Cb);
dy(1)=-(alpha*(3.5-y(2))/2/y(1))*To/T;
dy(2)=rate/Fao;
-------------------------------------

The ODE45 file is as follows:
-------------------------------------
[w,y]=ode45('ExampleGasPD_X_2013',[0 70],[1 0]);
figure; plot(w,y(:,1));
xlabel('w');
ylabel('P/Po');
figure; plot(w,y(:,2));
xlabel('W');
ylabel('X');

T=T1/T2;
figure; plot(y(:,2),T);
xlabel('w');
ylabel('T');
--------------------------------------
 
Status
Not open for further replies.

Part and Inventory Search

Sponsor

Back
Top