How to extract value at a given time in ode 45

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.

Answer 1

[Mischa Kim]

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.

%user parameters
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);
figure(1)
plot(t,y(:,2),t,y(:,4));
xlim(tspan);
ylabel('Population (n)');
xlabel('Time (days)');
legend('Infected','Location','SouthEast');

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);
end

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
end

[Shubham]

Thanks for helping Sir.

Answer 2

[KALYAN ACHARJYA]

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: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=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

t=0:5;
x=t.^2;
plot(t,x,'*');
hold on;

Next You may have use polyfit function to fit the data points, please refer MATLAB Doc

p=polyfit(t,x,7);
%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

t1=linspace(0,5);
x1=polyval(p,t1);
plot(t1,x1);

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.

t_req=2.3;
[near_val,idx]=min(abs(t1-t_req))

Results:

near_val =
    0.0232
idx =
    47

Once you get the idx, get the corresponding x values from the fitting data points, hence x1 value at idx index position

x1(idx)

Result:

ans =
    5.4222

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).

[Shubham]

Thanks for helping Sir.