Сложные значения возвращаются через некоторое время при решении ODE с SciPy - PullRequest
Новая
       104

Сложные значения возвращаются через некоторое время при решении ODE с SciPy

0 голосов
/ 27 января 2019

Я моделирую движение двух планет, вращающихся вокруг звезды, в присутствии диска (мелких частиц пыли / газа), распространяющегося по некоторому радиусу от звезды.Планеты становятся все ближе и ближе к звезде в центре.Рядом со звездой находится (пустая) полость с радиусом rcav в коде.Способ моделирования этого состоит в том, чтобы ввести расстояние (от rcav до R), на котором определенные значения (временные рамки демпфирования ta1, ta2, te1, te2) увеличиваются до очень больших значений (я использовал 10e + 12), то есть временные масштабы (ось y) изменяется как this на расстояниях a1 и a2 (ось x) (на графике R = 5 и rcav = 2).

Приведенный ниже код хорошо работает довремя t достигает значения около 2.82945e + 06.На этом этапе значения ta1 и dn1dt (и, возможно, другие) возвращают комплексные значения, и решатель ODE останавливается.

Как видите, я работаю с очень грязными ODE.Я также извиняюсь за грязный код, но я никогда не использовал Python до этого.Первоначально, ODE не смог бы преодолеть время t = 2e + 06.Я думаю, что это было связано с большим количеством я работал с.Итак, я узнал о пакете mpmath, который, казалось, работал довольно хорошо.Я также использовал тождество x = exp (log (x)), чтобы «упростить» некоторые уравнения для решателя.

Я также попытался решить эту проблему в R. Там же я столкнулся с проблемой, когда решатель будетзакончить шаги («превышены maxsteps») до окончания. 

import numpy as np
from scipy.integrate import odeint
import mpmath as mp
mp.dps = 10; mp.pretty = True

#specify constants
q = 1
m1 = 1*6*10**(24)
m2 = 1.6*6*10**(24)
mearth = 6*10**(24)
msun = 2*10**30
h = 0.02
md = 3*10**(-5)*2*10**30
mstar = 2*10**30
a = 3.5/5.26
rcav = 0.2
R = 0.22
tinf = 10**12

f1 = -1.3913023689045727
f2 = 0.5392267278220553

#define ODEs
def g(z, t):
    n1 = z[0]
    n2 = z[1]
    e1 = z[2]
    e2 = z[3]
    vp1 = z[4]
    vp2 = z[5]
    sig = z[6]
    a1 = z[7]
    a2 = z[8]
    d = z[9]
    phi1 = sig - vp1
    phi2 = sig - vp2

    ta1R = 27*(1 + (e1/(1.3*h))**5)*(1 - (e1/(1.1*h))**4)**(-1)*(h/0.02)**2*(msun/md)*(mearth/m1)*R*(msun/mstar)**0.5
    ta2R = 27*(1 + (e2/(1.3*h))**5)*(1 - (e2/(1.1*h))**4)**(-1)*(h/0.02)**2*(msun/md)*(mearth/m2)*R*(msun/mstar)**0.5
    te1R = (3*10**(-2))*(1 + 0.25*(e1/h)**3)*(h/0.02)**4*(msun/md)*(mearth/m1)*R*(msun/mstar)**0.5
    te2R = (3*10**(-2))*(1 + 0.25*(e2/h)**3)*(h/0.02)**4*(msun/md)*(mearth/m2)*R*(msun/mstar)**0.5
    if (a1 > R and a2 > R):
      ta1 = 27*(1 + (e1/(1.3*h))**5)*(1 - (e1/(1.1*h))**4)**(-1)*(h/0.02)**2*(msun/md)*(mearth/m1)*a1*(msun/mstar)**0.5
      ta2 = 27*(1 + (e2/(1.3*h))**5)*(1 - (e2/(1.1*h))**4)**(-1)*(h/0.02)**2*(msun/md)*(mearth/m2)*a2*(msun/mstar)**0.5
      te1 = (3*10**(-2))*(1 + 0.25*(e1/h)**3)*(h/0.02)**4*(msun/md)*(mearth/m1)*a1*(msun/mstar)**0.5
      te2 = (3*10**(-2))*(1 + 0.25*(e2/h)**3)*(h/0.02)**4*(msun/md)*(mearth/m2)*a2*(msun/mstar)**0.5
    elif ((R > a1) and (a1 > rcav) and (a2 > R)):
      ta1 = mp.exp(mp.log(tinf) + mp.log(ta1R/tinf)*((mp.mpf(a1) - rcav)/(R - rcav)))
      ta2 = 27*(1 + (e2/(1.3*h))**5)*(1 - (e2/(1.1*h))**4)**(-1)*(h/0.02)**2*(msun/md)*(mearth/m2)*a2*(msun/mstar)**0.5
      te1 = mp.exp(mp.log(tinf) + mp.log(te1R/tinf)*((mp.mpf(a1) - rcav)/(R - rcav)))
      te2 = (3*10**(-2))*(1 + 0.25*(e2/h)**3)*(h/0.02)**4*(msun/md)*(mearth/m2)*a2*(msun/mstar)**0.5
    elif ((R > a1) and (a1 > rcav) and (R > a2) and (a2 > rcav)):
      ta1 = mp.exp(mp.log(tinf) + mp.log(ta1R/tinf)*((a1 - rcav)/(R - rcav)))
      ta2 = mp.exp(mp.log(tinf) + mp.log(ta2R/tinf)*((a2 - rcav)/(R - rcav)))
      te1 = mp.exp(mp.log(tinf) + mp.log(te1R/tinf)*((a1 - rcav)/(R - rcav)))
      te2 = mp.exp(mp.log(tinf) + mp.log(te2R/tinf)*((a2 - rcav)/(R - rcav)))
    else:
      ta1 = tinf
      ta2 = tinf
      te1 = tinf
      te2 = tinf

    dn1dt = -3*q*mp.mpf(n1)**2*(a*m2/mstar)*(mp.mpf(e1)*f1*mp.sin(phi1) + e2*f2*mp.sin(phi2)) + (3*mp.mpf(n1))/(2*ta1) + (3*mp.mpf(n1)*e1**2)/(te1)
    dn2dt = 3*(q+1)*n2**2*(m1/mstar)*(e1*f1*np.sin(phi1) + e2*f2*np.sin(phi2)) + (3*n2)/(2*ta2) + (3*n2*e2**2)/(te2)
    de1dt = -n1*(a*m2/mstar)*(f1*np.sin(phi1)) - e1/te1
    de2dt = -n2*(m1/mstar)*(f2*np.sin(phi2)) - e2/te2
    dvp1dt = n1*(a*m2/mstar)*(1/e1)*f1*np.cos(phi1)
    dvp2dt = n2*(m1/mstar)*(1/e2)*f2*np.cos(phi2)
    dsigdt = (q+1)*n2 - q*n1
    da1dt = -2/3*mp.mpf(a1)*mp.mpf(dn1dt.real)/mp.mpf(n1)
    da2dt = -2/3*a2*dn2dt/n2
    dddt = 3*(q+1)*(n2/n1)*(e1*f1*np.sin(phi1) + e2*f2*np.sin(phi2))*((q+1)*n2*(m1/mstar) + q*n1*(a*m2/mstar)) + (n2/n1)*(3*(q+1))*(1/(2*ta2) + e2**2/te2 - 1/(2*ta1) - e1**2/te1)
    dzdt = [dn1dt, dn2dt, de1dt, de2dt, dvp1dt, dvp2dt, dsigdt, da1dt, da2dt, dddt]

    print(t,ta1,ta2,a1,a2,dn1dt) #to see where the error occurs
    print("--------")
    return dzdt

#inital condition
z0 = [2*np.pi*(3.5)**(-3/2), q/(q+1)*2*np.pi*(3.5)**(-3/2), 0.0001, 0.00011, np.pi/3, np.pi/5, np.pi/7, 3.5, 5.26, -0.0001]

# time points
t = np.linspace(0,4*10**6,num = 5000)

# solve ODE
sol = odeint(g, z0, t, atol = 1e-10, rtol = 1e-10)

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

Заранее благодарен за помощь.

...