Błąd Odd SciPy ODE Integracja

Question

Wdrażam bardzo prosty model podatny na zakażenie - odzyskany, z ciągłą populacją dla nieużytecznego projektu pobocznego - zwykle dość banalnym zadaniem. Ale używam PysCeS lub SciPy, które wykorzystują lsoda jako podstawowy solver. Dzieje się tak tylko w przypadku określonych wartości parametru i nie jestem pewien dlaczego. Kod używam jest następująco:

import numpy as np 
from pylab import * 
import scipy.integrate as spi 

#Parameter Values 
S0 = 99. 
I0 = 1. 
R0 = 0. 
PopIn= (S0, I0, R0) 
beta= 0.50  
gamma=1/10. 
mu = 1/25550. 
t_end = 15000. 
t_start = 1. 
t_step = 1. 
t_interval = np.arange(t_start, t_end, t_step) 

#Solving the differential equation. Solves over t for initial conditions PopIn 
def eq_system(PopIn,t): 
    '''Defining SIR System of Equations''' 
    #Creating an array of equations 
    Eqs= np.zeros((3)) 
    Eqs[0]= -beta * (PopIn[0]*PopIn[1]/(PopIn[0]+PopIn[1]+PopIn[2])) - mu*PopIn[0] + mu*(PopIn[0]+PopIn[1]+PopIn[2]) 
    Eqs[1]= (beta * (PopIn[0]*PopIn[1]/(PopIn[0]+PopIn[1]+PopIn[2])) - gamma*PopIn[1] - mu*PopIn[1]) 
    Eqs[2]= gamma*PopIn[1] - mu*PopIn[2] 
    return Eqs 

SIR = spi.odeint(eq_system, PopIn, t_interval) 

To daje następujący błąd:

lsoda-- at current t (=r1), mxstep (=i1) steps 
     taken on this call before reaching tout  
     In above message, I1 =  500 
     In above message, R1 = 0.7818108252072E+04 
Excess work done on this call (perhaps wrong Dfun type). 
Run with full_output = 1 to get quantitative information. 

Normalnie, kiedy wystąpi problem takiego, coś śmiertelnie złego układu równań I utworzonej , ale nie widzę w tym nic złego. Dziwnie, działa również, jeśli zmienisz mu na coś podobnego do 1/15550. W przypadku było coś nie tak z systemem, ja również realizowane w modelu R następująco:

require(deSolve) 

sir.model <- function (t, x, params) { 
    S <- x[1] 
    I <- x[2] 
    R <- x[3] 
    with (
    as.list(params), 
{ 
    dS <- -beta*S*I/(S+I+R) - mu*S + mu*(S+I+R) 
    dI <- beta*S*I/(S+I+R) - gamma*I - mu*I 
    dR <- gamma*I - mu*R 
    res <- c(dS,dI,dR) 
    list(res) 
} 
) 
} 

times <- seq(0,15000,by=1) 
params <- c(
beta <- 0.50, 
gamma <- 1/10, 
mu <- 1/25550 
) 

xstart <- c(S = 99, I = 1, R= 0) 

out <- as.data.frame(lsoda(xstart,times,sir.model,params)) 

Ten również wykorzystuje lsoda, ale wydaje się, że dzieje się bez żadnych przeszkód. Czy ktoś może zobaczyć, co się dzieje w kodzie Pythona?

Answer 1

Myślę, że dla parametrów, które wybrałeś masz problemy z stiffness - ze względu na niestabilność numeryczną, rozmiar kroku solver staje się bardzo mały w regionach, w których nachylenie krzywej rozwiązania jest rzeczywiście całkiem płytki. Solver firmy Fortran lsoda, który jest opakowany przez scipy.integrate.odeint, próbuje adaptacyjnie przełączać się między metodami dopasowanymi do "sztywnych" i "nie-sztywnych" systemów, ale w tym przypadku wydaje się, że nie przechodzi na sztywne metody.

Bardzo grubsza można po prostu znacznie zwiększyć maksymalny dozwolony kroki i solver dostanie tam w końcu:

SIR = spi.odeint(eq_system, PopIn, t_interval,mxstep=5000000) 

Lepszym rozwiązaniem jest użycie obiektowego ODE solver scipy.integrate.ode, który pozwala jednoznacznie określić, czy użycie sztywnych lub non-sztywnych metod:

import numpy as np 
from pylab import * 
import scipy.integrate as spi 

def run(): 
    #Parameter Values 
    S0 = 99. 
    I0 = 1. 
    R0 = 0. 
    PopIn= (S0, I0, R0) 
    beta= 0.50  
    gamma=1/10. 
    mu = 1/25550. 
    t_end = 15000. 
    t_start = 1. 
    t_step = 1. 
    t_interval = np.arange(t_start, t_end, t_step) 

    #Solving the differential equation. Solves over t for initial conditions PopIn 
    def eq_system(t,PopIn): 
     '''Defining SIR System of Equations''' 
     #Creating an array of equations 
     Eqs= np.zeros((3)) 
     Eqs[0]= -beta * (PopIn[0]*PopIn[1]/(PopIn[0]+PopIn[1]+PopIn[2])) - mu*PopIn[0] + mu*(PopIn[0]+PopIn[1]+PopIn[2]) 
     Eqs[1]= (beta * (PopIn[0]*PopIn[1]/(PopIn[0]+PopIn[1]+PopIn[2])) - gamma*PopIn[1] - mu*PopIn[1]) 
     Eqs[2]= gamma*PopIn[1] - mu*PopIn[2] 
     return Eqs 

    ode = spi.ode(eq_system) 

    # BDF method suited to stiff systems of ODEs 
    ode.set_integrator('vode',nsteps=500,method='bdf') 
    ode.set_initial_value(PopIn,t_start) 

    ts = [] 
    ys = [] 

    while ode.successful() and ode.t < t_end: 
     ode.integrate(ode.t + t_step) 
     ts.append(ode.t) 
     ys.append(ode.y) 

    t = np.vstack(ts) 
    s,i,r = np.vstack(ys).T 

    fig,ax = subplots(1,1) 
    ax.hold(True) 
    ax.plot(t,s,label='Susceptible') 
    ax.plot(t,i,label='Infected') 
    ax.plot(t,r,label='Recovered') 
    ax.set_xlim(t_start,t_end) 
    ax.set_ylim(0,100) 
    ax.set_xlabel('Time') 
    ax.set_ylabel('Percent') 
    ax.legend(loc=0,fancybox=True) 

    return t,s,i,r,fig,ax 

wyjściowa:

Answer 2

Zarażona populacja PopIn[1] rozpada się do zera. Najwyraźniej (normalna) niedokładność liczbowa prowadzi do uzyskania wartości ujemnej (około -3.549e-12) w pobliżu t = 322,9. W końcu roztwór rozwinie się w pobliżu t = 7818.093, z PopIn[0] idącym w kierunku + nieskończoności i PopIn[1] zbliżającym się do nieskończoności.

Edytuj: Usunąłem wcześniejszą sugestię dotyczącą "szybkiej poprawki". To było wątpliwe włamanie.