时间戳

首页 > 解决方案 > 如何解决常微分方程组解的数值不稳定性

numerical-methods - 如何解决常微分方程组解的数值不稳定性

问题描述

我一直在尝试获得以下常微分方程组的数值解：

倾斜午餐中身体在空气中的运动方程：（
显然 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.00001if(fabs(theta)<0.00001) theta = -0.00001theta

标签： numerical-methodsdifferential-equationsapproximationrunge-kutta

解决方案

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

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

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

推荐阅读