24 views (last 30 days)

Show older comments

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

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

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

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!