How to extract value at a given time in ode 45
112 views (last 30 days)
I have solved an ode numerically for tspan [0 5] so an array is created for values of x at different t. The question is hoe can I extract value of x near say 2 and print it. I am asking this because I need to plot x(2.3) against a parameter.
Mischa Kim on 18 Jan 2021
Hi Shubham, use an events functions, see the example below. With events functions you can identify zero crossings; in your case (when t = 2.3) the zero crossing you would want to detect is for the expression: t - 2.3. In the command where ode45 is called below MATLAB returns te and ye the time of the event and the value of the function at this time.
N = 45742000; %total population
I0 = 1; %initial infected population
r0 = 12.9; %reproductive value
R0 = 0;%initial recovered population
i_period = 9; %duration of infectious period
tspan = [1, 50]; %length of time to run simulation over
%rate constant calculations
mu = 1 / i_period; %recovery rate
beta = r0 * mu / N; %infection rate
S0 = N - I0 - R0; %initial susceptible population
N0 = I0 + R0 + S0; %total population
%feeding parameters into function
pars = [beta, mu, N, r0];
y0 = [45742000 1 0 0];
%using the ode45 function to perform intergration
options = odeset('Events',@events);
[t,y,te,ye,ie] = ode45(@sir_rhs, tspan, y0, options, pars);
function f = sir_rhs(~,y,pars)
f = zeros(4,1);
f(1) = -pars(1)*y(1)*y(2);
f(2) = pars(1)*y(1)*y(2) - pars(2)*y(2);
f(3) = pars(2) * y(2);
f(4) = y(2);
function [value,isterminal,direction] = events(t,~,~)
value = t - 20; % This would be your 2.3 (t - 2.3)
isterminal = 0; % Do not stop the integration when zero is detected
direction = 0; % Detect all zeros
More Answers (1)
KALYAN ACHARJYA on 18 Jan 2021
Edited: KALYAN ACHARJYA on 18 Jan 2021
Hope I have understood your question, if not, please do not hesitate to correct me.
Lets understand with an example, say t as array with range 0 to 5
t =[0 1 2 3 4 5]
And you have evaluated the x, which is function of t, lets say x=t^2
x=[0 1 4 9 16 25]
Now, you want to get the x value at t = 2.3, right? Lets plot the original data
%7>Use polyfit to fit a 7th-degree polynomial to the points.
% Please see the other Polynomial order also
Create the data range as fiting datapoints
See the results, original data points and fiting datapoints
Next, you want to get the x value the same as t = 2.3, x(2.3), In MATLAB or any other languages, you have to carefully deals with floating points numbers. Most cases 1.0 is not equal to 1. Refer the links, for more details Link1, Link2
Next the task is, get the indix position (idx) of t1, where t1 value is nearest to 2.3, hence I have used absolate (irrespective of sign) and get the minimum, please check the near_val, which is still not zero in subtraction cases, which menas the t1 datapoints is not exactly as 2.3, see other poly order and practice more.
Once you get the idx, get the corresponding x values from the fitting data points, hence x1 value at idx index position
Hope this helps, I know, the answer is quite long, but it is a very important issue in the number of cases, so that it can help others too (novices).