Simulating ODEs with array inputs

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;

Answer 1

[Stephane Dauvillier]

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 ?

[Muhammad Jibran Ejaz]

Thanx mate. It solved the problem like Thanos' snap.

[Jan]

A linear interpolation of the parameters is still not smooth. While ode45 might compute a final value, the intergrator is driven apart from its specifications. For a numerical simulation of a scientific problem, this is not reliable.

Answer 2

[Jan]

Of course the integrator must crash. Check the used dimensions: What is the size of

dydt = (eta*P+(T_o-y)/R)/C;

when you define the parameters as row vectors? The implicit expanding creates matrices, and as the error message tells you, the output of the function to be integrated must be a column vector.

You did not mention, what you want to achieve. I guess, that T_o is a parameter, which cheanges its value every second. Then first_order() is a non-smooth function and not compatible with Matlab's ODE integrators. Remember, that ODE45 is designed for smooth functions only, because otherwise the stepsize controller is driven out of its specifications. The calculated integral can be dominated by rounding or discreatization erros, such that it is as useful as a random value. Don't do this in scientific work.

If a parameter is changed in different intervals, you have to run the integration in these intervals and use the final value of one interval as initial value of the next one:

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 = 1:length(T_o);
...
allT = 1;
allY = y0;
for it = 1:length(T_o)
  aT_o   = T_o(it);
  t      = it;
  [T, Y] = ode45(@(t, y)first_order(t, y, P, T_o(t)), [t, t+1], y0);
  allT   = [allT; T(2:end)];
  allY   = [allY; Y(2:end, :)];
  y0     = Y(end, :);
end

[Muhammad Jibran Ejaz]

Thanx. This works for me as well.