Как решить численную неустойчивость к решению системы обыкновенных дифференциальных уравнений - PullRequest
Новая
       80

Как решить численную неустойчивость к решению системы обыкновенных дифференциальных уравнений

0 голосов
/ 11 декабря 2018

Я пытался получить численное решение для следующей системы обыкновенных дифференциальных уравнений:

Уравнения для движения тела по воздуху в наклонном обеде:
(очевидно, LaTeX не делаетработа над переполнением стека) 

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

по алгоритму Рунге-Кутта-Фельберга, но в середине вычисления мне нужно вычислить тета, который вычисляется как 

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 попадает в ловушку этого отрицательного значения, есть ли у кого-нибудь указания на то, что нужно сделать, чтобы устранить эту числовую нестабильность?

1 Ответ

0 голосов
/ 11 декабря 2018

Плохая идея использовать обратные косинусные или синусоидальные функции для определения угла точки.Чтобы получить 

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

чтобы вам не понадобился какой-либо угол или попытка ослабить связь компонентов скорости, уменьшив ее до угла вектора скорости.Теперь вы можете подключиться к выбранному вами методу интеграции, применяемому как метод для векторных систем.

...