how to change the coefficient of system of ODE with respect to time ?

조회 수: 5 (최근 30일)
i recently found a matlab code in https://in.mathworks.com/. i have enclosed that code in this section. in this code, i want to change the coefficient value (transm and recov) with respect to time as given in excel file . Is it possible to change this code with respect to time ?
% Just some start parameters and coefficients
Istart = 0.01;
Sstart = 0.99;
Rstart = 0;
transm = 3.2;
recov = 0.23;
maxT = 20;
% define S',I',R'
% S' = - transm * S * I
% I' = transm * S I - recov * I
% R' = recov * I
% S is y(1)
% I is y(2)
% R is y(3)
%so here fun = @(t,y)[S'; I'; R'];
fun = @(t,y)[-transm*y(1)*y(2); (transm*y(1)*y(2))-(recov*y(2)); recov*y(2)];
% Provide the starting parameters
Y0 = [Sstart; Istart; Rstart;];
% Define the range of t
tspan = [0 maxT];
% Magic happens and matrix Y contains S,I,R
[T,Y] = ode45(fun,tspan,Y0);
% Plot plot
figure
plot(T,abs(Y(:,1)),'b-') % S
hold on
plot(T,abs(Y(:,2)),'r-') % I
plot(T,abs(Y(:,3)),'g-') % R

채택된 답변

Star Strider
Star Strider 2020년 9월 10일
Interpolate from the Excel file, and choose a stiff solver:
% Just some start parameters and coefficients
Istart = 0.01;
Sstart = 0.99;
Rstart = 0;
maxT = 20;
T1 = readtable('data.xls');
interptransm = @(t) interp1(T1.t, T1.transm, t, 'linear','extrap'); % Interpolate
interprecov = @(t) interp1(T1.t, T1.recov, t, 'linear','extrap'); % Interpolate
%so here fun = @(t,y)[S'; I'; R'];
fun = @(t,y)[-interptransm(t)*y(1)*y(2); (interptransm(t)*y(1)*y(2))-(interprecov(t)*y(2)); interprecov(t)*y(2)];
% Provide the starting parameters
Y0 = [Sstart; Istart; Rstart;];
% Define the range of t
tspan = [0 maxT];
% Magic happens and matrix Y contains S,I,R
[T,Y] = ode15s(fun,tspan,Y0);
% Plot plot
figure
semilogy(T,abs(Y(:,1)),'b-') % S
hold on
plot(T,abs(Y(:,2)),'r-') % I
plot(T,abs(Y(:,3)),'g-') % R
hold off
grid
legend('S', 'I', 'R', 'Location', 'N')
producing:
.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

태그

Community Treasure Hunt

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

Start Hunting!

Translated by