How to extract value at a given time in ode 45
12 views (last 30 days)
Show older comments
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.
0 Comments
Accepted Answer
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.
%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
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: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;
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/491310/image.png)
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).
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!