如何求解常微分方程组解的数值不稳定性

Question

我一直试图得到以下常微分方程组的数值解：

body 在倾斜的午餐中通过空气运动的方程式：
（显然 LaTeX 对堆栈溢出不起作用）

u'= -F(u, theta, t)*cos(theta)
v'= -F(v, theta, t)*sin(theta)-mg

通过Runge-Kutta-Fehlberg算法，但在计算的中间我必须计算theta，即通过

计算

arccos(u/sqrt(u^2+v^2)) or arcsin(v/sqrt(u^2+v^2)),

但最终 theta 变得太小了，我需要它来求解函数 F( v, theta, t) 并找到值 V = sqrt(v^2 + u^2) 我使用 V = (v/sin(theta))，但是作为 theta 变小，sin(theta) 也变小，我从给定的向前迭代中得到一个数值错误 -1.IND00，这可能是因为 theta 太小了，我试图使 theta 从像 0.00001 这样的小正角度到像 -0.00001 (if(fabs(theta)<0.00001) theta = -0.00001) 这样的小负角度，但似乎 theta 陷入了这个负值，有没有人有关于如何解决这种数值不稳定性的指示？

Answer 1

使用反余弦或正弦函数来确定点的角度是个坏主意。得到

theta = arg ( u + i*v)

使用

theta = atan2(v,u).

这仍然存在跳到负半轴的问题，即v=0, u<0。这可以通过使 theta 成为第三个动态变量来解决，因此

 theta' = Im( (u'+i*v')/(u+i*v) ) = (u*v' - u'*v) / (u^2+v^2)

但实际上，具有空气摩擦力的自由落体方程最容易实现为

def friction(vx, vy):
    v = hypot(vx, vy)
    return k*v

def freefall_ode(t, u):
    rx, ry, vx, vy = u
    f=friction(vx, vy)
    ax = -f*vx
    ay = -f*vy - g
    return array([ vx, vy, ax, ay ])

这样您就不需要任何角度或尝试通过将其减小到速度矢量的角度来削弱速度分量的耦合。您现在可以将其插入您选择的积分方法中，作为矢量值系统的一种方法应用。

如何求解常微分方程组解的数值不稳定性

How to solve the numerical instability to a solution of a system of ordinary differential equations

numerical-methods

approximation

differential-equations

runge-kutta