Would like to graph an ode45 answer with if statment

2 次查看（过去 30 天）
Rodrigo Blas
Rodrigo Blas 2020-4-19
评论： Walter Roberson 2020-4-20
function diffy2=Ljallo_f_q2(~,X)
y=X(1);T=X(2);
%%data
Q=0.06555; %%m^3/s
ctot=0.0203; %%knol.m^3
alpha=26900; %%m^2/m^3
V=6*10^-4; %%m^3
Mhat=30; %%kg/kmol
cpg=1070; %%J/kg*k
delhrxn=-2.84*10^8; %%j/kmol
ep=0.68; %%unitless
cps=1000; %%j/kg*k
rhog=ctot*Mhat; %%kg/m^3
rhos=2500; %%kg/m^3
if t<1
yin=.01;Tin=600;
if (1<t)&&(t<2)
yin=.03;
Tin=700;
else
yin=.03;
Tin=700;
end
end
zspan=[0 2];
ic=[yin Tin]
k1=6.70*10^10*exp(-12556/T);
K1=65.5*exp(961/T);
r=(0.05*k1*y)/(T(1+K1*y)^2);
diffy2(1)=(Q*ctot*(yin-y)-alpha*V*r)/(ep*V*ctot);
diffy2(2)=Q*ctot*Mhat*cpg*(Tin-T)+alpha*V*(-delhrxn)*r/(V*(ep*rhog*cpg+(1-ep)*rhos*cps));
diffy2=diffy2';
end
>> [t,y]=ode45(@Ljallo_f_q2,zspan,ic)
Undefined function or variable 't'.
Error in Ljallo_f_q2 (line 16)
if t<1
Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets
args{1} to yp0.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode,
tspan, y0, options, varargin);
>>
I would like to graph this ODE with the if statemetn. The graph goes from 0 to 2 on the x axis(t) with the if conditions between 0 and 2. The initial conditions change based on the x-axis(t).I understand that my t is not defined but I dont know how to arrange it so that matlab reads it properly.

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

回答（1 个）

Walter Roberson
Walter Roberson 2020-4-19
Change
function diffy2=Ljallo_f_q2(~,X)
to
function diffy2=Ljallo_f_q2(t,X)
Also you will need to move
zspan=[0 2];
ic=[yin Tin];
into your script that runs ode45 instead of having them inside your Ljallo_f_q2
Do not be surprised if ode45 gives you completely wrong answers, or complains about not being able to find a solution at time 1. The mathematics used by the ode*() functions require that your equations are continuous to two derivatives, but your equations are not continuous at t==1 or t==2.
What you should be doing is running ode45 from time 0 to 1*(1-eps), then stopping it and using the output as the initial conditions to running ode45 from time 1 to 2*(1-eps), and then stopping it and using the output as the initial conditions to run time 2 exactly.
  2 个评论
Rodrigo Blas
Rodrigo Blas 2020-4-19
function diffy2=Ljallo_f_q2(t,X)
y=X(1);T=X(2);
%%data
Q=0.06555; %%m^3/s
ctot=0.0203; %%knol.m^3
alpha=26900; %%m^2/m^3
V=6*10^-4; %%m^3
Mhat=30; %%kg/kmol
cpg=1070; %%J/kg*k
delhrxn=-2.84*10^8; %%j/kmol
ep=0.68; %%unitless
cps=1000; %%j/kg*k
rhog=ctot*Mhat; %%kg/m^3
rhos=2500; %%kg/m^3
yin=.01;
Tin=600;
k1=6.70*10^10*exp(-12556/T);
K1=65.5*exp(961/T);
r=(0.05*k1*y)/(T*(1+K1*y)^2);
diffy2(1)=(Q*ctot*(yin-y)-alpha*V*r)/(ep*V*ctot);
diffy2(2)=Q*ctot*Mhat*cpg*(Tin-T)+alpha*V*(-delhrxn)*r/(V*(ep*rhog*cpg+(1-ep)*rhos*cps));
diffy2=diffy2';
end
yin=.01;
Tin=600;
zspan=[0 1];
ic=[yin Tin];
[t,y]=ode45(@Ljallo_f_q2,zspan,ic);
plot(t,y(:,1));
xlabel('time');
ylabel('Mol Fraction CO');
Thank you for fixing that for me. I set it from 0 to 1 I'm getting a straight hoorizontal line now. Which is incorrect.
How would I fix this?
Walter Roberson
Walter Roberson 2020-4-20
You are not getting a straight horizontal line. You have again made the mistake of not accounting for scale.
Give the command
xlim([0.005 1]); ylim auto
and you will see that the line is not straight at all!

请先登录，再进行评论。

类别

在 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