使用fsolve和Euler 3BDF的SIR模型 [英] SIR model using fsolve and Euler 3BDF

问题描述

您好，我被要求使用MATLAB中的fsolve命令来求解SIR模型，而Euler 3向后指向.我对如何继续感到很困惑，请帮忙.这是我到目前为止所拥有的.我为3BDF方案创建了一个函数，但是我不确定如何进行fsolve和求解非线性ODE系统. SIR模型显示为，而3BDF方案则表示为

Hi i've been asked to solve SIR model using fsolve command in MATLAB, and Euler 3 point backward. I'm really confused on how to proceed, please help. This is what i have so far. I created a function for 3BDF scheme but i'm not sure how to proceed with fsolve and solve the system of nonlinear ODEs. The SIR model is shown as and 3BDF scheme is formulated as

clc
clear all 
gamma=1/7;
beta=1/3;
ode1= @(R,S,I) -(beta*I*S)/(S+I+R);
ode2= @(R,S,I) (beta*I*S)/(S+I+R)-I*gamma;
ode3= @(I) gamma*I;
f(t,[S,I,R]) = [-(beta*I*S)/(S+I+R); (beta*I*S)/(S+I+R)-I*gamma; gamma*I];
R0=0;
I0=10;
S0=8e6;

odes={ode1;ode2;ode3}
fun = @root2d;
x0 = [0,0];
x = fsolve(fun,x0)



function [xs,yb] = ThreePointBDF(f,x0, xmax, h, y0)
% This function should return the numerical solution of y at x = xmax.
% (It should not return the entire time history of y.)
% TO BE COMPLETED


xs=x0:h:xmax;
y=zeros(1,length(xs));
y(1)=y0;
yb(1)=y0+f(x0,y0)*h;


for i=1:length(xs)-1

R =R0;


y1(i+1,:) = fsolve(@(u) u-2*h/3*f(t(i+1),u) - R, y1(i-1,:)+2*h*F(i,:))


S = S0;
y2(i+1,:) = fsolve(@(u) u-2*h/3*f(t(i+1),u) - S, y2(i-1,:)+2*h*F(i,:))


I= I0;
y3(i+1,:) = fsolve(@(u) u-2*h/3*f(t(i+1),u) - I, y3(i-1,:)+2*h*F(i,:))



end


end

推荐答案

您有一个隐式方程

y(i+1) - 2*h/3*f(t(i+1),y(i+1)) = R = (4*y(i) - y(i-1))/3

，其中右侧项R在该步骤中是恒定的.

where the right-side term R is constant in the step.

请注意，这是针对矢量值系统y'(t)=f(t,y(t))其中

f(t,[S,I,R]) = [-(beta*I*S)/(S+I+R); (beta*I*S)/(S+I+R)-I*gamma; gamma*I];

以某种方式.

要解决此问题

R = (4*y(i,:) - y(i-1,:))/3
y(i+1,:) = fsolve(@(u) u-2*h/3*f(t(i+1),u) - R, y(i-1,:)+2*h*F(i,:))

其中中点步长用于获得2阶近似作为初始猜测.如有必要，添加求解器选项以提高容错能力.一个人也只能保留一小段函数值，然后必须注意该短数组中的位置与时间索引的对应关系.

where a midpoint step is used to get an order 2 approximation as initial guess. Add solver options for error tolerances if necessary. One can also keep only a short array of function values, then one has to be careful for the correspondence of the position in the short array to the time index.

使用所有这些和参考解决方案，由高阶标准求解器生成以下组件误差图

Using all that and a reference solution by a higher order standard solver produces the following error graphs for the components

在其中可以看到，恒定的第一步的一阶误差导致一阶全局误差，而在第一步中使用Euler方法产生二阶误差会导致明显的二阶全局误差.

where one can see that the first order error of the constant first step results in a first order global error, while with a second order error in the first step using the Euler method results in a clear second order global error.

实施一般方法

from scipy.optimize import fsolve

def BDF2(f,t,y0,y1):
    N = len(t)-1;
    y = (N+1)*[np.asarray(y0)];
    y[1] = y1;
    for i in range(1,N):
        t1, R = t[i+1], (4*y[i]-y[i-1])/3
        y[i+1] = fsolve(lambda u: u-2*h/3*f(t1,u)-R, y[i-1]+2*h*f(t[i],y[i]), xtol=1e-3*h**3)
    return np.vstack(y)

设置要求解的模型

gamma=1/7;
beta=1/3;
print beta, gamma
y0 = np.array([8e6, 10, 0])
P = sum(y0); y0 = y0/P
def f(t,y): S,I,R = y; trns = beta*S*I/(S+I+R); recv=gamma*I; return np.array([-trns, trns-recv, recv])

为两个初始化变量计算参考解决方案和方法解决方案

Compute a reference solution and method solutions for the two initialization variants

from scipy.integrate import odeint

tg = np.linspace(0,120,25*128)
yg = odeint(f,y0,tg,atol=1e-12, rtol=1e-14, tfirst=True)

M = 16; # 8,4
t = tg[::M];
h = t[1]-t[0];
y1 = BDF2(f,t,y0,y0)
e1 = y1-yg[::M]
y2 = BDF2(f,t,y0,y0+h*f(0,y0))
e2 = y2-yg[::M]

按照上面的方法绘制误差，但可以将其嵌入到plot命令中，原则上可以通过首先计算解决方案列表来进行分离

Plot the errors, computation as above, but embedded in the plot commands, could be separated in principle by first computing a list of solutions

fig,ax = plt.subplots(3,2,figsize=(12,6))
for M in [16, 8, 4]:
    t = tg[::M];
    h = t[1]-t[0];
    y = BDF2(f,t,y0,y0)
    e = (y-yg[::M])
    for k in range(3): ax[k,0].plot(t,e[:,k],'-o', ms=1, lw=0.5, label = "h=%.3f"%h)
    y = BDF2(f,t,y0,y0+h*f(0,y0))
    e = (y-yg[::M])
    for k in range(3): ax[k,1].plot(t,e[:,k],'-o', ms=1, lw=0.5, label = "h=%.3f"%h)
for k in range(3): 
    for j in range(2): ax[k,j].set_ylabel(["$e_S$","$e_I$","$e_R$"][k]); ax[k,j].legend(); ax[k,j].grid()
ax[0,0].set_title("Errors: first step constant");
ax[0,1].set_title("Errors: first step Euler")

这篇关于使用fsolve和Euler 3BDF的SIR模型的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！