- The ode solver should be called as oriented object fully, i think that help
- The ode function must have linked the parameter to gamma and beta
- The ode funcion must return the values of the slopes (it was fixed at dydt = [0;0;0])
SENSITIVITY ANALYSIS OF A SYSTEM OF AN ODE USING NORMALIZED SENSITIVITY FUNCTION
22 ビュー (過去 30 日間)
古いコメントを表示
I want to perform the sensitivity analysis on the parameters of an ODE SIR model using normalized sensitivity function. I used the following code below. When it was run, I got this error: "undefined functionnor variable 'p', error in epidemic1 line8, F=ode45(@epidemic,ic,p)". How can I resolve it to get the plots?
My interest is to get the figure (3) i.e. The plot of nomalized sensitivity functions against time
function epidemic1()
clc
beta=0.4; gamma=0.04;
p1 = [beta gamma];
ic = [995;5;0];
F = ode45(@epidemic,ic,p);
sol1 = solve(F,0,80)
figure(1)
plot(sol1.Time,sol1.Solution,'-o')
legend('S','I','R')
title(['SIR Populations over Time','$\beta=0.4$, $\gamma=0.04$'],'Interpreter','latex')
xlabel('Time','Interpreter','latex')
ylabel('Population','Interpreter','latex')
Figure(2)
p2 = [0.2 0.1];
F.Parameters = p2;
F.Sensitivity = odeSensitivity;
sol2 = solve(F,0,80)
plot(sol2.Time,sol2.Solution,'-o')
legend('S','I','R')
title(['SIR Populations over Time','\beta=0.2$, $\gamma=0.1$'],'Interpreter','latex')
xlabel('Time','Interpreter','latex')
ylabel('Population','Interpreter','latex')
U11 = squeeze(sol2.Sensitivity(1,1,:))'.*(p2(1)./sol2.Solution(1,:));
U12 = squeeze(sol2.Sensitivity(1,2,:))'.*(p2(2)./sol2.Solution(1,:));
U21 = squeeze(sol2.Sensitivity(2,1,:))'.*(p2(1)./sol2.Solution(2,:));
U22 = squeeze(sol2.Sensitivity(2,2,:))'.*(p2(2)./sol2.Solution(2,:));
U31 = squeeze(sol2.Sensitivity(3,1,:))'.*(p2(1)./sol2.Solution(3,:));
U32 = squeeze(sol2.Sensitivity(3,2,:))'.*(p2(2)./sol2.Solution(3,:));
Figure(3)
t = tiledlayout(3,2);
title(t,['Parameter Sensitivity by Equation','$\beta=0.2$, $\gamma=0.1$'],'Interpreter','latex')
xlabel(t,'Time','Interpreter','latex')
ylabel(t,'\% Change in Eqn','Interpreter','latex')
nexttile
plot(sol2.Time,U11)
title('$\beta$, Eqn. 1','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U12)
title('$\gamma$, Eqn. 1','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U21)
title('$\beta$, Eqn. 2','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U22)
title('$\gamma$, Eqn. 2','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U31)
title('$\beta$, Eqn. 3','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U32)
title('$\gamma$, Eqn. 3','Interpreter','latex')
ylim([-5 5])
function dydt = epidemic(t,y,p)
dydt = [0;0;0];
S = y(1);
I = y(2);
R = y(3);
N = S + I + R;
dSdt = -beta*(I*S)/N;
dIdt = beta*(I*S)/N - gamma*I;
dRdt = gamma*I;
end
end
2 件のコメント
Simon Diaz
2024 年 11 月 1 日 1:33
I found some issues. Here are the main :
Also figures must be called with lowercase (figure not Figure)
%% Corrected version epidemic1()
clc;close all;clear all
beta=0.4; gamma=0.04;
p1 = [beta gamma];
ic = [995;5;0];
fun_ode = @(t,x,p)epidemic(t,x,p);
problem_ode = ode( ODEFcn = fun_ode ); % Define ODE function
problem_ode.InitialValue = ic; % Inivial value x(t=t0)
problem_ode.Parameters = p1;
sol1 = solve(problem_ode,0,80)
figure(1)
plot(sol1.Time,sol1.Solution,'-o')
legend('S','I','R')
title(['SIR Populations over Time','$\beta=0.4$, $\gamma=0.04$'],'Interpreter','latex')
xlabel('Time','Interpreter','latex')
ylabel('Population','Interpreter','latex')
figure(2)
p2 = [0.2 0.1];
problem_ode.Parameters = p2;
problem_ode.Sensitivity = odeSensitivity;
sol2 = solve(problem_ode,0,80)
plot(sol2.Time,sol2.Solution,'-o')
legend('S','I','R')
title(['SIR Populations over Time','\beta=0.2$, $\gamma=0.1$'],'Interpreter','latex')
xlabel('Time','Interpreter','latex')
ylabel('Population','Interpreter','latex')
U11 = squeeze(sol2.Sensitivity(1,1,:))'.*(p2(1)./sol2.Solution(1,:));
U12 = squeeze(sol2.Sensitivity(1,2,:))'.*(p2(2)./sol2.Solution(1,:));
U21 = squeeze(sol2.Sensitivity(2,1,:))'.*(p2(1)./sol2.Solution(2,:));
U22 = squeeze(sol2.Sensitivity(2,2,:))'.*(p2(2)./sol2.Solution(2,:));
U31 = squeeze(sol2.Sensitivity(3,1,:))'.*(p2(1)./sol2.Solution(3,:));
U32 = squeeze(sol2.Sensitivity(3,2,:))'.*(p2(2)./sol2.Solution(3,:));
figure(3)
t = tiledlayout(3,2);
title(t,['Parameter Sensitivity by Equation','$\beta=0.2$, $\gamma=0.1$'],'Interpreter','latex')
xlabel(t,'Time','Interpreter','latex')
ylabel(t,'\% Change in Eqn','Interpreter','latex')
nexttile
plot(sol2.Time,U11)
title('$\beta$, Eqn. 1','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U12)
title('$\gamma$, Eqn. 1','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U21)
title('$\beta$, Eqn. 2','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U22)
title('$\gamma$, Eqn. 2','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U31)
title('$\beta$, Eqn. 3','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U32)
title('$\gamma$, Eqn. 3','Interpreter','latex')
ylim([-5 5])
function dydt = epidemic(t,y,p)
S = y(1);
I = y(2);
R = y(3);
N = S + I + R;
beta = p(1);
gamma = p(2);
dSdt = -beta*(I*S)/N;
dIdt = beta*(I*S)/N - gamma*I;
dRdt = gamma*I;
dydt = [dSdt;dIdt;dRdt];
end
回答 (1 件)
Star Strider
2024 年 6 月 9 日
I am not certain that this does everything you want, however it now has the virtue of running without error —
epidemic1
function epidemic1()
clc
beta=0.4; gamma=0.04;
p1 = [beta gamma];
ic = [995;5;0];
F = ode;
F.ODEFcn = @(t,y)epidemic(t,y,p1);
F.InitialValue = ic;
F.SelectedSolver
sol1 = solve(F,0,80)
figure(1)
plot(sol1.Time,sol1.Solution,'-o')
legend('S','I','R')
title(["SIR Populations over Time","$\beta=0.4$, $\gamma=0.04$"],'Interpreter','latex')
xlabel('Time','Interpreter','latex')
ylabel('Population','Interpreter','latex')
figure(2)
p2 = [0.2 0.1];
F.Parameters = p2;
F.Sensitivity = odeSensitivity;
sol2 = solve(F,0,80)
plot(sol2.Time,sol2.Solution,'-o')
legend('S','I','R')
title(['SIR Populations over Time','$\beta=0.2$, $\gamma=0.1$'],'Interpreter','latex')
xlabel('Time','Interpreter','latex')
ylabel('Population','Interpreter','latex')
U11 = squeeze(sol2.Sensitivity(1,1,:))'.*(p2(1)./sol2.Solution(1,:));
U12 = squeeze(sol2.Sensitivity(1,2,:))'.*(p2(2)./sol2.Solution(1,:));
U21 = squeeze(sol2.Sensitivity(2,1,:))'.*(p2(1)./sol2.Solution(2,:));
U22 = squeeze(sol2.Sensitivity(2,2,:))'.*(p2(2)./sol2.Solution(2,:));
U31 = squeeze(sol2.Sensitivity(3,1,:))'.*(p2(1)./sol2.Solution(3,:));
U32 = squeeze(sol2.Sensitivity(3,2,:))'.*(p2(2)./sol2.Solution(3,:));
figure(3)
t = tiledlayout(3,2);
title(t,["Parameter Sensitivity by Equation","$\beta=0.2$, $\gamma=0.1$"],'Interpreter','latex')
xlabel(t,'Time','Interpreter','latex')
ylabel(t,'\% Change in Eqn','Interpreter','latex')
nexttile
plot(sol2.Time,U11)
title('$\beta$, Eqn. 1','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U12)
title('$\gamma$, Eqn. 1','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U21)
title('$\beta$, Eqn. 2','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U22)
title('$\gamma$, Eqn. 2','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U31)
title('$\beta$, Eqn. 3','Interpreter','latex')
ylim([-5 5])
nexttile
plot(sol2.Time,U32)
title('$\gamma$, Eqn. 3','Interpreter','latex')
ylim([-5 5])
function dydt = epidemic(t,y,p)
dydt = [0;0;0];
S = y(1);
I = y(2);
R = y(3);
N = S + I + R;
dSdt = -beta*(I*S)/N;
dIdt = beta*(I*S)/N - gamma*I;
dRdt = gamma*I;
end
end
.
5 件のコメント
Star Strider
2024 年 6 月 9 日
@SAHEED AJAO — It would be best to upgrade. However, you can run it here (although not on MATLAB Online, since it reflects your version and installed Toolboxes).
As @Sam Chak mentioned, you can calculate it manually. If you have the Symbolic Math Toolbox, using the jacobian function and then the matlabFunction function makes that easier.
Sam Chak
2024 年 6 月 9 日
The formula can be found here:
参考
カテゴリ
Help Center および File Exchange で Ordinary Differential Equations についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!