python - 使用欧拉方法求解微分方程组

Question

我正在尝试用欧拉方法求解常微分方程组，但是当我尝试打印速度时，我得到了

RuntimeWarning: overflow encountered in double_scalars

我得到的不是打印数字而是 nan(不是数字)。我认为问题可能出在定义加速度时，但我不确定，如果有人能帮助我，我将不胜感激。

from numpy import *
from math import pi,exp

d=0.0005*10**-6
a=[]

while d<1.0*10**-6 :
    d=d*2
    a.append(d)

D=array(a)

def a_particula (D,x, v, t):
    cc=((1.00+((2.00*lam)/D))*(1.257+(0.400*exp((-1.10*D)/(2.00*lam)))))
    return (g-((densaire*g)/densparticula)-((mu*18.0*v)/(cc*densparticula*  (D**2.00))))

def euler (acel,D, x, v, tv, n=15):
    nv, xv, vv = tv.size, zeros_like(tv), zeros_like(tv)
    xv[0], vv[0] = x, v
    for k in range(1, nv):
        t, Dt = tv[k-1], (tv[k]-tv[k-1])/float(n)
        for i in range(n):
            a = acel(D,x, v, t)
            t, v = t+Dt, v+a*Dt
            x = x+v*Dt
        xv[k], vv[k] = x, v
    return (xv, vv)

g=(9.80)    
densaire= 1.225         
lam=0.70*10**-6         
densparticula=1000.00
mu=(1.785*10**-5)       
tv = linspace(0, 5, 50)
x, v = 0, 0 #initial conditions


for j in range(len(D)):

    xx, vv = euler(a_particula, D[j], x, v, tv)

    print(D[j],xx,vv)

Answer 1

如果您在问题中包含完整的警告消息，将来会有所帮助 - 它将包含出现问题的行:

tmp/untitled.py:15: RuntimeWarning: overflow encountered in double_scalars
  return (g-((densaire*g)/densparticula)-((mu*18.0*v)/(cc*densparticula*  (D**2.00))))

Overflow当变量的大小超过可以表示的最大值时发生。在本例中，double_scalars 指的是一个 64 位 float ，其最大值为:

print(np.finfo(float).max)
# 1.79769313486e+308

所以表达式中有一个标量值:

(g-((densaire*g)/densparticula)-((mu*18.0*v)/(cc*densparticula*  (D**2.00))))

超过 ~1.79e308。要找出是哪一个，您可以使用 np.errstate在发生这种情况时引发 FloatingPointError，然后捕获它并启动 Python debugger :

    ... 
    with errstate(over='raise'):
        try:
            ret = (g-((densaire*g)/densparticula)-((mu*18.0*v)/(cc*densparticula*  (D**2.00))))
        except FloatingPointError:
            import pdb
            pdb.set_trace()
        return ret
    ...

然后您可以在调试器中检查此表达式各个部分的值。溢出似乎发生在:

(mu*18.0*v)/(cc*densparticula*  (D**2.00))

第一次出现警告时，(cc*densparticula* (D**2.00) 计算为 2.3210168586496022e-12，而 (mu*18.0*v)计算结果为 -9.9984582297025182e+299。

基本上你是用一个非常大的数字除以一个非常小的数字，结果的大小超过了可以表示的最大值。这可能是您的数学问题，也可能是您对函数的输入未合理缩放。