Solving simultaneous ODEs using ode45
6 views (last 30 days)
Show older comments
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 Comments
Accepted Answer
Torsten
on 22 Jan 2016
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 Comments
Torsten
on 25 Jan 2016
I don't like giving the same name to variables which have different meanings.
I think that the main problem with your code was that your variable lists in your call to f33trial were inconsistent.
Best wishes
Torsten.
More Answers (0)
See Also
Categories
Find more on Performance and Memory in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!