Метод Рунге-Кутты четвертого порядка. Интеграция в обратном направлении - PullRequest
Новая
       41

Метод Рунге-Кутты четвертого порядка. Интеграция в обратном направлении

3 голосов
/ 14 октября 2019

Я использую метод Рунге-Кутты четвертого порядка для численного решения обычного уравнения движения фонового скалярного поля в искривленном пространстве-времени с квартичным потенциалом: 

$\phi^{''}=-3\left(1+\frac{H^{'}}{3H}\right)\phi^{'}-\lambda\phi^3/H^2$,

$'$ denoting the derivative w.r.t. the e-folds number $\textrm{d}N=H\textrm{d}t$ and, from the Friedmann equation:

$H^2=\frac{\lambda \phi^4}{4}\frac{1}{3M_{Pl}^2-(1/2)\phi^{'2}}$;

$H^{'}=-\frac{1}{2M_{Pl}^2}H\phi^{'2}$.

Проблема возникает при интегрировании в обратном направлении с использованиемВ качестве начальных условий конечные значения я получил после интегрирования вперед. Результат взорвать, не сопоставляя значения, полученные ранее, при интеграции вперед. Я просто не понимаю, в чем проблема, поскольку и уравнение, и код не являются неизвестными. Во-первых, я интегрировал от 0 до 64 электронных складок. Тогда я просто поменяю направление интеграции.

Я также прилагаю код: 

def rk4trial(f,v0,t0,tf,n,V):  
    t=np.linspace(t0,tf,n)
    h=t[1]-t[0]
    v=np.array((n+1)*[v0])
    for j in range(n):  
        V.append(v[j])
        V[j]=copy.deepcopy(V[j])
        k1=f(v[j],t[j])*h
        k2=f(v[j]+(1/2)*k1,t[j]+(1/2)*h)*h
        k3=f(v[j]+(1/2)*k2,t[j]+(1/2)*h)*h
        k4=f(v[j]+k3,t[j]+h)*h
        v[j+1]=v[j]+(k1+2*k2+2*k3+k4)/6
    return V, t, h


def Fdet(v,t):
    phi, sigma = v
    H=(((lamb/4)*phi**4)/(3*mpl**2-(1/2)*sigma**2))**(1/2)
    HH=-((1/2)*sigma**2)*(1/mpl**2)
    return np.array([sigma,-3*(1+HH/3)*sigma-lamb*phi**3/(H**2)])

PS: Этот вопрос также был размещен здесь: https://scicomp.stackexchange.com/questions/33583/runge-kutta-fourth-order-method-integrating-backwards,, где уравнения показаны подробно.

РЕДАКТИРОВАТЬ: ненужные части в коде удалены.

РЕДАКТИРОВАТЬ: В ответ на @LutzL, я присоединяю графики \ phi / M_ {Pl} и \ phi ^ {'} после интегрирования вперед (сплошные линии) и назад (пунктирные линии), выполнивчто она сказала. Как видите, есть резкое отклонение от результатов прямой интеграции, которое я не могу объяснить.

|\phi^{'}| vs N

\phi/M_{Pl} vs N

1 Ответ

1 голос
/ 15 октября 2019

Я бы изменил метод RK4 до минимально необходимого. Нет необходимости иметь массив v, частично дублирующий содержимое массива V, поэтому 

def rk4trial(f,v0,t0,tf,n,V):  
    t=np.linspace(t0,tf,n)
    h=t[1]-t[0]
    v=v0
    for j in range(n):  
        V.append(v)
        k1=f(v,t[j])*h
        k2=f(v+0.5*k1,t[j]+0.5*h)*h
        k3=f(v+0.5*k2,t[j]+0.5*h)*h
        k4=f(v+k3,t[j]+h)*h
        v=v+(k1+2*k2+2*k3+k4)/6
    return V, t, h

Здесь нет проблем с копированием, поскольку v создается заново на каждом шаге,так что все объекты, добавляемые в возвращаемый массив, являются отдельными.

Обратная интеграция должна быть такой же простой, как прямая интеграция, 

V1, t1, h1 = rk4trial(Fdet,v0,t0,tf,n,[])
V2, t2, h2 = rk4trial(Fdet,V1[-1],tf,t0,n,[])

и V2[k] должны совпадать с V1[-k-1] в пределах ошибки метода. Большие различия следует ожидать только в жестких ODE.

...