Solving simultaneous ODEs using ode45

Question

XXC 2016-1-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/264720-solving-simultaneous-odes-using-ode45

评论： Torsten 2016-1-25

I am trying to solve 4 simultaneous ODEs using ode45. One of these coupled ODE is rather complicated, with a heaviside time dependance. I am not very familiar with Matlab so I'm not sure what went wrong when Error messages came up.

The function is:

if true
  % code
end
function output = f33trial(V,F,rho,Cp,Tc0,a,b,T,Fc,Fcs,dFc,T0s,t)
% solving 4 simultaneous ODEs in part 3.3
output=zeros(4,1);
Fc=Fcs+dFc*heaviside(t-5);
output(1)=F/0.7/V*(T0s-T(1))-Fc/0.1/V*(T(1)-Tc0)*(1-exp(-a/rho/Cp*Fc^(b-1)));
output(2)=F/0.1/V*(T(1)-T(2));
output(3)=F/0.1/V*(T(2)-T(3));
output(4)=F/0.1/V*(T(3)-T(4));
end

Here is the script:

if true
  % code
end
% Defining constants
a=5.9e5; % J/(min^0.5dCm^1.5)
b=0.5;
Cp=4187; % J/kgdC
rho=1000; % kg/m^3
F=0.085; % m^3/min
Tc0=25; % dC
V=2.1; % m^3
T0s=150; % dC
Tsps=85; % dC
dFc=0.05; % m^3/min
Fcs=0.51775; % m^3/min obtained from Script32
tspan=[0,10];
Tstart=[T0s,T0s,T0s,T0s];
options=odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
[t,T]=ode45(@(t,T) f33trial(T,t),tspan,Tstart,options);
plot(t,T)

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

请先登录，再进行评论。

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

Answer 1

Torsten 2016-1-22

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/264720-solving-simultaneous-odes-using-ode45#answer_207036

在 MATLAB Online 中打开

Try this:

function main
% Defining constants
a=5.9e5; % J/(min^0.5dCm^1.5)
b=0.5;
Cp=4187; % J/kgdC
rho=1000; % kg/m^3
F=0.085; % m^3/min
Tc0=25; % dC
V=2.1; % m^3
T0s=150; % dC
Tsps=85; % dC
dFc=0.05; % m^3/min
Fcs=0.51775; % m^3/min obtained from Script32
tspan=[0,10];
Tstart=[T0s,T0s,T0s,T0s];
options=odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
[T,Y]=ode45(@(t,y) f33trial(t,y,a,b,Cp,rho,F,Tc0,V,T0s,Tsps,dFc,Fcs),tspan,Tstart,options);
plot(T,Y)
function output = f33trial(t,T,a,b,Cp,rho,F,Tc0,V,T0s,Tsps,dFc,Fcs)
output=zeros(4,1);
Fc=Fcs+dFc*heaviside(t-5);
output(1)=F/0.7/V*(T0s-T(1))-Fc/0.1/V*(T(1)-Tc0)*(1-exp(-a/rho/Cp*Fc^(b-1)));
output(2)=F/0.1/V*(T(1)-T(2));
output(3)=F/0.1/V*(T(2)-T(3));
output(4)=F/0.1/V*(T(3)-T(4));
end

Best wishes

Torsten.

2 个评论
显示无隐藏无

XXC 2016-1-22

Thanks :) It works now but I'm not sure why the original didn't. Also, what is the variable Y (and y) that you put in and why does it make a difference?