MATLAB: Solving simultaneous ODEs using ode45

ode45

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)

Best Answer

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.

Related Solutions

MATLAB: F=@(t,y) 4exp(0.8t)-0.5*y; tspan = [0 4]; y0 = 2; h = 1; [t,y] = ode_euler(f,tspan,y0,h)

tspan = [0 4];
y0 = 2;
[t1,y1] = ode45(@(t,y) 4*exp(0.8*t)-0.5*y,tspan,y0,h);
figure('name','Reda','color','w');
plot(t1,y1, 'r-','LineWidth',2);grid on;
xlabel('Time');ylabel('Solution y(Time)');
title('And you instead of a blanket substitute solver Euler');

MATLAB: How to write the code for open luup nd closed luup

The answer for open luup is :-

clc
close all
clear all
h = 0.5; %in minutes
run_count = 201; %time count
%FIXED PARAMETERS
cp = 1; %cal/gm C

rho = 1e6; %gm/cubic meters

cpc = 1; %cal/gm C
rhoc = 1e6; %gm/cubic meters
a = 1.41e5;
b = 0.5;
Tcs = 25;
%VARYING PARAMETERS
v = 3.5; %cubic meters
fs = 0.085; %cubic meters/min

fcs = 0.44; %cubic meters/min
Ths= 150;
alp = fs/v;
beta1 = a*fcs^(b+1)/(rho*cp*v);
beta2 = a*fcs^b/(2*rhoc*cpc);
beta = beta1/(fcs + beta2);
Ts = (alp*Ths + beta*Tcs)/(alp + beta);
Tc = Tcs; %constant
Th = Ths*ones(run_count,1); %disturbance

f = fs*ones(run_count,1); %disturbance
fc = fcs*ones(run_count,1); %manipulated variable
t = zeros(run_count,1);
T = Ts*ones(run_count,1);
noise_std = 0.07;
Tm = zeros(run_count,1);
Tm(1) = T(1) + noise_std*randn;
%introduction of disturbance
n = input('Specify the variable for introducing step change. \n Enter 1 for coolant flow rate. \n Enter 2 for inlet stream temperature.\n');
step_time = 21;
if n == 1
    delta_fc = 0.1*fcs; % introducing +10% change in flow rate
    fc(step_time:end) = delta_fc + fc(step_time:end);
else
    delta_t = +10; %introducing a +10 deg. C change in inlet temperature
    Th(step_time:end) = Th(step_time:end) + delta_t;
end
for k = 1:run_count-1
    alp_ = f(k)/v;
    beta1_ = a*fc(k)^(b+1)/(rho*cp*v);
    beta2_ = a*fc(k)^b/(2*rhoc*cpc);
    beta_ = beta1_/(fc(k) + beta2_);    
    T(k+1) = T(k) + h*[alp_*(Th(k) - T(k)) - beta_*(T(k) - Tc)];    
    Tm(k+1) = T(k+1) + noise_std*randn;
    t(k+1)=k*0.5;
end
set(0,'DefaultLineLineWidth',2)
set(0,'DefaultaxesLineWidth',2)
set(0,'DefaultaxesFontSize',16)
set(0,'DefaultTextFontSize',16)
set(0,'DefaultaxesFontName','arial')
figure(1), plot(t, T,'b', t, Tm, 'g',t, Ts*ones(run_count,1),'r'), grid;
xlabel('Time (min)')
ylabel('STHE Temperature (deg C)')
title('Controlled Output (T)')
legend('True Temperature','Measured Temperature','Initial S.S. Temperature')
figure(2), stairs(t, fc), grid;
xlabel('Time (min)')
ylabel('Coolant Flow (cubic meters/min)')
title('Manipulated Input (F_c)')
figure(3), stairs(t, Th), grid;
xlabel('Time (min)')
ylabel('Inlet Stream Temperature (deg C)')
title('Disturbance (T_h)')
diff = [mean(Tm(end-10:end))-mean(Tm(1:10))];
delta = 0.63*diff*ones(run_count,1);
figure(4), plot(t, Tm - Ts,'b' ,t, delta,'r'), grid;
xlabel('Time (min)')
ylabel('Deviation STHE Temperature (deg C)')
title('Deviation Controlled Output (T)')
legend('Measured Output','63.2% of Change in Output')
figure(5), stairs(t, fc-fcs), grid;
xlabel('Time (min)')
ylabel('Deviation Coolant Flow (cubic meters/min)')
title('Deviation Manipulated Input (Fc)')
figure(6), stairs(t, Th - Ths), grid;
xlabel('Time (min)')
ylabel('Deviation Inlet Temperature (deg C)')
title('Deviation Disturbance (T)')
clc
if n == 1
    K_mani = diff/delta_fc
else
    K_dist = diff/delta_t
end

Related Question

Solving coupled differential equations