Azzera filtri
Azzera filtri

SENSITIVITY ANALYSIS OF A SYSTEM OF AN ODE USING NORMALIZED SENSITIVITY FUNCTION

11 visualizzazioni (ultimi 30 giorni)
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

Risposte (1)

Star Strider
Star Strider il 9 Giu 2024
I am not certain that this does everything you want, however it now has the virtue of running without error —
epidemic1
ans =
SolverID enumeration ode45
sol1 =
ODEResults with properties: Time: [0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 62 64 66 68 70 72 74 76 78 80] Solution: [3x41 double]
sol2 =
ODEResults with properties: Time: [0 0.0027 27.4608 80] Solution: [3x4 double] Sensitivity: [3x2x4 double]
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 Commenti
Star Strider
Star Strider il 9 Giu 2024
@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.

Accedi per commentare.

Prodotti


Release

R2013a

Community Treasure Hunt

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

Start Hunting!

Translated by