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

Question

SAHEED AJAO 2024-6-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2126826-sensitivity-analysis-of-a-system-of-an-ode-using-normalized-sensitivity-function

评论： Sam Chak 2024-6-9

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2024-6-9

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2126826-sensitivity-analysis-of-a-system-of-an-ode-using-normalized-sensitivity-function#answer_1469221

在 MATLAB Online 中打开

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

.

5 个评论
显示 3更早的评论隐藏 3更早的评论

SAHEED AJAO 2024-6-9

在 MATLAB Online 中打开

Thanks, I am grateful. I modified the code and ran it , I used Matlab R2013a to run it. I got this error: "Undefined function or variable 'ode'.

Error in epidemic11 (line 7)

F = ode; "

The modified code is given below

function epidemic11()
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 yprime = epidemic(~,y,~)
   yprime = [0;0;0]; 
    S = y(1);
    I = y(2);
    R = y(3);
    N = S + I + R;
    yprime(1) = -beta*(I*S)/N;
    yprime(2) = beta*(I*S)/N - gamma*I;
    yprime(3) = gamma*I;
end
end

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.