基于fsolve和Euler 3BDF的SIR模型

Question

嗨，我被要求使用MATLAB中的fsolve命令来解决SIR模型，并向后返回Euler 3点。我真的不知道该怎么做，请帮助我。这就是我到目前为止所拥有的。我为3BDF格式创建了一个函数，但我不确定如何继续求解非线性常微分方程组。SIR模型如下所示

并且3BDF方案被表示为

​

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

Answer 1

你有一个隐式方程

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

其中，右侧项G在对fsolve的调用中是常量，即在求解隐式步长方程期间。

请注意，这适用于向量值系统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];

要解决此问题，请编写

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

其中，中点步长用于获得作为初始猜测的2阶近似，F(i,:)=f(t(i),y(i,:))。根据需要添加误差公差的求解器选项，您希望隐式方程中的误差小于步骤的截断误差O(h^3)。也可以只保留函数值的短数组，然后必须注意短数组中的位置与时间索引的对应关系。

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

​

​

其中可以看到，常数第一步的一阶误差导致一阶全局误差，而使用欧拉方法的第一步中的二阶误差导致明显的二阶全局误差。

在一般情况下实现方法

from scipy.optimize import fsolve

def BDF2(f,t,y0,y1):
    N, h = len(t)-1, t[1]-t[0];
    y = (N+1)*[np.asarray(y0)];
    y[1] = y1;
    for i in range(1,N):
        t1, G = t[i+1], (4*y[i]-y[i-1])/3
        y[i+1] = fsolve(lambda u: u-2*h/3*f(t1,u)-G, 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])

计算两个初始化变量的参考解和方法解

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命令中，原则上可以通过首先计算一组解决方案来分离

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")