Simulating ODEs with array inputs
2 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Muhammad Jibran Ejaz
il 29 Apr 2019
Commentato: Jan
il 29 Apr 2019
I have an equation that I want to simulate:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/217159/image.png)
The P is calculated by the following equation:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/217160/image.png)
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;
0 Commenti
Risposta accettata
Stephane Dauvillier
il 29 Apr 2019
Modificato: Jan
il 29 Apr 2019
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 ?
2 Commenti
Jan
il 29 Apr 2019
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.
Più risposte (1)
Jan
il 29 Apr 2019
Modificato: Jan
il 29 Apr 2019
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
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!