使用 fsolve 和 Euler 3BDF 的 SIR 模型

Question

您好，我被要求在 MATLAB 中使用 fsolve 命令求解 SIR 模型，并且 Euler 3 指向倒退。我真的很困惑如何继续，请帮助。这就是我到目前为止所拥有的。我为 3BDF 方案创建了一个函数，但我不确定如何继续使用 fsolve 并求解非线性 ODE 系统。 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)) where

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]

绘制错误，计算如上，但嵌入在绘图命令中，原则上可以通过首先计算解决方案列表来分离

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