用欧拉方法求解微分方程组

Question

我试图用欧拉方法求解一个常微分方程组，但当我试图打印速度时，我得到

RuntimeWarning: overflow encountered in double_scalars

而不是打印数字，我得到nan（不是数字）。我想问题可能出在定义加速度的时候，但我不确定，如果有人能帮我，我会非常感激的。在

^{pr2}$

Answer 1

你的系统减少到

dv/dt = a = K - L*v

其中K大约10和{}之间，乍一看1e+5到{}。使用的实际系数证实：

^{pr2}$

---

速度的欧拉步长是

v[j+1]=v[j]+(K-L*v[j])*dt =(1-L*dt)*v[j] + K*dt

对于任何类似于预期的摩擦效应的东西，即速度降到K/L，我们需要abs(1-L*dt)<1，如果可能{}，也就是dt < 1/L。这意味着dt < 1e-10。在

为了能够使用更大的时间步长，您需要使用stiff微分方程的方法，这意味着隐式方法。最简单的方法是隐式欧拉法、中点法和梯形法。在

由于线性关系，中点法和梯形法相当于同一公式

v[j+1] = v[j] + dt * ( K - L*(v[j]+v[j+1])/2 )

或者

v[j+1] = ( (1-L*dt/2)*v[j] + K*dt ) / (1+L*dt/2)

---

当然，最简单的方法是精确地集成ODE

(-L*v')/(K-L*v)=-L  =>  K-L*v(t)=C*exp(-L*t), C=K-L*v(0)

v(t)=K/L + exp(-L*t)*(v(0)-K/L)

集成到

x(t)=x(0)+K/L*t+(1-exp(-L*t))/L*(v(0)-K/L).

---

或者有可能你在把物理定律转换成公式时犯了一个错误，因此常数的大小都是错误的。在

Answer 2

将来，如果您在问题中包含完整的警告消息，这将很有帮助—它将包含发生问题的行：

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位浮点，其最大值为：

^{pr2}$

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

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

超过了~1.79e308。要找出哪一个，可以使用^{}在发生这种情况时引发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。在

基本上你把一个很大的数除以一个很小的数，结果的大小超过了可以表示的最大值。这可能是您的数学问题，也可能是您对函数的输入没有合理地缩放。在