How to extract value at a given time in ode 45

조회 수: 29 (최근 30일)
Shubham 2021년 1월 18일
편집: KALYAN ACHARJYA 2021년 1월 18일
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 2021년 1월 18일
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
댓글 수: 1이전 댓글 -1개 표시이전 댓글 -1개 숨기기
Shubham 2021년 1월 18일
Thanks for helping Sir.

댓글을 달려면 로그인하십시오.

추가 답변 (1개)

KALYAN ACHARJYA 2021년 1월 18일
편집: KALYAN ACHARJYA 2021년 1월 18일
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).
댓글 수: 1이전 댓글 -1개 표시이전 댓글 -1개 숨기기
Shubham 2021년 1월 18일
Thanks for helping Sir.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Programming에 대해 자세히 알아보기

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by