MATLAB: Simulating ODEs with array inputs

MATLABodeode45

I have an equation that I want to simulate:
The P is calculated by the following equation:
The code I've tried is as follows:
clear all; close all; clc
T_o = [26 25 24 23 22 22 21 23 24 26 27 28 29 30 31 32 31 29 29 27 26 25 24 22];
T_in = 26;
eta = 2.6;
R = 5.56;
P_avg = (T_out - T_in)/(eta*R);
P_avg(T_out <= T_in) = 0;
P_avg(T_out >= T_in + 2) = 1490;
t = 1:length(T_o);
y0 = 28;
P = P_avg;
[T Y] = ode45(@(t, y)first_order(t, y, P, T_o), t, y0);
The function I've defined is as follows:
function dydt = first_order(t, y, P, T_o)
% function dydt = first_order(t, y, fP, fT_o)
% P = interp1(fP, t, t);
% T_o = interp1(fT_o, t, t);
R = 5.56;
eta = 2.6;
C = 0.18;
dydt = (eta*P+(T_o-y)/R)/C;
end
The command window gives the following error:
Error using odearguments (line 90)
@(T,Y)FIRST_ORDER(T,Y,P,T_O) must return a column vector.
Error in ode45 (line 113)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...
Error in simulate (line 25)
[T Y] = ode45(@(t, y)first_order(t, y, P, T_o), t, y0);
When I try the code with single values, it works fine. But as try to pass the whole array of T_o and P, it crashs and generates the above mentioned error. I have tried the solution in https://www.mathworks.com/help/matlab/ref/ode45.html#examples
as you can see in the commented parts but somehow I can't make it work. Any help will be appreciated;

Best Answer

Hi, I'm not sure you use ode45 correctly.
In your function dydt = first_order(t, y, P, T_o)
t is the current time and y the current value, P and T_o are parameters of your function.
Meaning that it won't only be equal to 1 or 2 or .... or length(T_o).
If I've understood correctly the parameters P and T_o depends on t so P and T_o should be function.
So in your main code, you may want to have this:
fcnTout = @(time) interp1(t,T_out,time);
fcnP = @(time) interp1(t,P_avg,time);
[T,Y] = ode45(@(t, y)first_order(t, y, fcnP, fcnTout), t, y0);
figure;plot(T,Y); % just to see the result
And in your first_order function:
dydt = (eta*P(t)+(T(t)_o-y)/R)/C;
Does it solve your problem ?
Related Question