Odd SciPy Lỗi tích hợp ODE

Question

Tôi đang triển khai mô hình Nhạy cảm-Bị nhiễm-Phục hồi rất đơn giản với dân số ổn định cho một dự án bên nhàn rỗi - thường là một nhiệm vụ khá tầm thường. Nhưng tôi đang chạy vào lỗi giải quyết bằng cách sử dụng hoặc PysCeS hoặc SciPy, cả hai đều sử dụng lsoda như là người giải quyết cơ bản của họ. Điều này chỉ xảy ra đối với các giá trị cụ thể của một tham số và tôi được cho là lý do tại sao. Mã Tôi đang sử dụng như sau:

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) 

này tạo ra các lỗi sau:

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. 

Thông thường khi tôi gặp phải một vấn đề như thế, có điều gì đó ở giai đoạn cuối xảy ra với hệ thống phương trình tôi thiết lập , nhưng cả hai tôi không thể nhìn thấy bất cứ điều gì sai trái với nó. Thật kỳ lạ, nó cũng hoạt động nếu bạn thay đổi mu thành một cái gì đó như 1/15550. Trong trường hợp đó là một cái gì đó xảy ra với hệ thống này, tôi cũng thực hiện các mô hình trong R như sau:

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

này cũng sử dụng lsoda, nhưng có vẻ là đi ra mà không có một xô. Bất cứ ai có thể thấy những gì đang xảy ra sai trong mã Python?

Answer 1

Tôi nghĩ rằng đối với các thông số bạn đã chọn bạn đang gặp phải sự cố với stiffness - do tính không ổn định về số, kích thước bước của trình giải quyết được đẩy trở thành rất nhỏ ở các khu vực mà độ dốc của đường cong giải pháp thực sự khá nông. Bộ giải Fortran lsoda, được bọc bởi scipy.integrate.odeint, cố gắng chuyển đổi thích nghi giữa các phương pháp phù hợp với hệ thống 'cứng' và 'không cứng', nhưng trong trường hợp này có vẻ như không chuyển sang phương pháp cứng.

Rất thô sơ, bạn có thể chỉ ồ ạt tăng các bước tối đa cho phép và giải sẽ đạt được điều đó cuối cùng:

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

Một lựa chọn tốt hơn là sử dụng ODE hướng đối tượng giải scipy.integrate.ode, cho phép bạn dứt khoát chọn xem có nên sử dụng phương pháp cứng hoặc không cứng:

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 

Output:

Answer 2

Dân số bị nhiễm PopIn[1] phân rã thành 0. Rõ ràng, (bình thường) số không chính xác dẫn đến PopIn[1] trở thành âm (xấp xỉ -3.549e-12) gần t = 322.9. Sau đó, cuối cùng các giải pháp thổi lên gần t = 7818.093, với PopIn[0] đi về phía + vô cùng và PopIn[1] đi về phía-vô cực.

Chỉnh sửa: Tôi đã xóa đề xuất trước đó của mình về "sửa lỗi nhanh". Đó là một hack đáng ngờ.