使用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)) = R = (4*y(i) - y(i-1))/3

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

注意，这是针对矢量值系统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阶近似作为初始猜测。如有必要，添加求解器选项以提高容错能力。一个人也只能保留一小段函数值，然后必须注意该短数组中的位置与时间索引的对应关系。]