Scipy的Odeint在初始时间步解后返回零

Question

由于某种原因，在试图对ODE系统进行数值求解时，它会变成全零。我认为这意味着不稳定的解决方案，但总的来说，我认为情况将返回无穷大或nan，而不是零。为什么会发生这种情况，我该如何解决？我在运行odeint时收到以下警告（尽管它仍在完成）：

Warning

python\python37\lib\site-packages\scipy\integrate\odepack.py:236: ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information. 


 warnings.warn(warning_msg, ODEintWarning)

代码

import scipy as sp
from scipy.integrate import odeint
import numpy as np

def ODEEq(rho, t, D, B, Nd):
    return (D*np.gradient(np.gradient(rho))-B*(Nd*rho+(rho**2)))
def Solver(delz, delt, L, T, G0, D, B, Nd, Sf, Sb, unitMaker=1e-9):
    xm1Ghost = G0[1] - (Sf*delz*G0[0]/D)
    xLp1Ghost = G0[-2] - (Sb*delz*G0[-1]/D)
    rho = G0
    np.insert(rho, 0, xm1Ghost)
    np.append(rho, xLp1Ghost)

    sol = odeint(ODEEq, rho, tspace, args=(D, B, Nd), full_output=1)
    return sol

L = 4000
timeRange = 2000
D = 1.639e18
B = 2e13
Nd = 0
delz = L/(10*L)
delt = timeRange/(10*timeRange)
G0 = 640*np.exp((-6.4e-2)*np.linspace(0,1000,L))
Sf = 0.64
Sb = 0.64

sol = Solver(delz, delt, L, timeRange, G0, D, B, Nd, Sf, Sb)

结果（sol）：

(array([[6.40000000e+02, 6.29838965e+02, 6.19839253e+02, ...,
    1.05982469e-25, 1.04299825e-25, 1.02643897e-25],
   [3.64489112e+01, 3.64037414e+01, 3.63580081e+01, ...,
    1.64712980e-25, 1.61424855e-25, 1.58213241e-25],
   [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
    0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
   ...,
   [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
    0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
   [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
    0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
   [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, ...,
    0.00000000e+00, 0.00000000e+00, 0.00000000e+00]]), {'hu': array([2.24035532e-005, 8.91443222e-312, 4.78098028e+004, ...,
   3.21768899e-021, 3.11632803e-021, 3.01816006e-021]), 'tcur': array([1.28018934e-003, 8.91443222e-312, 2.65545333e-021, ...,
   4.31714361e-049, 4.15062501e-049, 3.98753839e-049]), 'tolsf': array([ 0.00000000e+000,  8.91443222e-312, -4.85590765e+004, ...,
    5.32418663e-009,  5.23965656e-009,  5.15646853e-009]), 'tsw': array([2.54756216e-004, 8.91443222e-312, 4.83671586e-009, ...,
   5.52645557e-023, 5.40880326e-023, 5.29122707e-023]), 'nst': array([        500,         420,  2051404512, ...,  -901543794,
    -382763050, -1396204354], dtype=int32), 'nfe': array([     52989,        420, 2051404512, ...,  896695326, -858716839,
   -564525364], dtype=int32), 'nje': array([        13,        420,  604126608, ..., 1834441325, 1211124015,
   1163472949], dtype=int32), 'nqu': array([         5,        420,  414912656, ..., -441395200, 1046325818,
    744409088], dtype=int32), 'imxer': -1, 'lenrw': 16036022, 'leniw': 4020, 'mused': array([         2,        420,  414912656, ...,   67999744, 1022034837,
    973886464], dtype=int32), 'message': 'Excess work done on this call (perhaps wrong Dfun type).'})

Answer 1

您的系统似乎是具有边界条件u_t = D*u_xx - B*(Nd*u+u**2)和D*u_x(0)=Sf*u(0)/2的D*u_x(L)=Sb*u(L)/2，这是具有非线性局部动力学的扩散过程，该非线性局部动力学在u=0处具有稳定的平衡，在u=-Nd处具有不稳定的平衡。

• 由于每个实现都有自由的边界，或者至少没有指定的边界，因此，系统在任何地方都收敛到其稳定的平衡并不是什么难解的。

• 您的设置令人费解，空间域的长度是多少，其细分是什么？您设置了delz = 0.1，但是np.linspace(0,1000,L)中的间距是1000/3999，大约是0.25

• 您不将空间分割应用于梯度计算，您将需要将delz作为参数传递，或者随后除以其平方。

• 您需要应用边界条件，分别鬼值的计算，以进行ODE函数的每次评估。

为了更深入地讨论如何将PDE转换为线法系统，请在scientific computing SE中发布到目前为止的方程式和您的想法，因为此处的数学内容不合时宜。据我所知，没有编程错误，只有这个可能构造错误的ODE系统的数理特征。