SENSITIVITY ANALYSIS OF A SYSTEM OF AN ODE USING NORMALIZED SENSITIVITY FUNCTION
23 次查看(过去 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
0 个评论
回答(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!