Модель SIR с использованием fsolve и Euler 3BDF

Question

Привет, меня попросили решить модель SIR с помощью команды fsolve в MATLAB, и Эйлер на 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 в качестве начального предположения. При необходимости добавьте параметры решателя для допусков ошибок. Можно также хранить только короткий массив значений функций, тогда нужно следить за соответствием позиции в коротком массиве временному индексу.

Используя все это и эталонное решение более высокого порядка Standard Solver создает следующие графики ошибок для компонентов

, где видно, что ошибка первого порядка постоянного первого шага приводит к первому порядку глобальная ошибка, в то время как ошибка второго порядка на первом шаге с использованием метода Эйлера приводит к явной глобальной ошибке второго порядка.

---

Реализуйте метод в общих чертах

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

Вычисление эталонного решения и методов решения для двух вариантов инициализации

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