CODEWITHSUNDEEP
matlabnumerical-methodsodedifferential-equations
Mv2156
Mv2156

Reputation: 27

SIR model using fsolve and Euler 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 enter image description here and 3BDF scheme is formulated as enter image description here

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

Upvotes: 1

Views: 629

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

You have an implicit equation

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

where the right-side term G is constant in the call to fsolve, that is, during the solution of the implicit step equation.

Note that this is for the vector valued system 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];

To solve this write

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,:))

where a midpoint step is used to get an order 2 approximation as initial guess, F(i,:)=f(t(i),y(i,:)). Add solver options for error tolerances as necessary, you want the error in the implicit equation smaller than the truncation error O(h^3) of the step. 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

enter image description here

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.

Implement the method in general terms

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)

Set up the model to be solved

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

Upvotes: 1

Related Questions