Runge-Kutta fourth order method. Integrating backwards

Question

I am using a Runge-Kutta fourth order method to solve numerically the usual equation of motion of a background scalar field in curved spacetime with a quartic potential:

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

$'$ denoting the derivative w.r.t. the e-folds number $	extrm{d}N=H	extrm{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}$.

The problem comes when integrating backwards using as the initial conditions the final values I got after integrating forward. The outcome blow up without matching the values obtained before, when integrating forward. I simply do not understand where the problem is as both the equation and the code are not unknown whatsoever. Firstly, I integrated from 0 to 64 e-folds. Then I simply reverse the integration direction.

I attach the code too:

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: This question has been posted here too:https://scicomp.stackexchange.com/questions/33583/runge-kutta-fourth-order-method-integrating-backwards, where equations are shown in detail.

EDIT: Unnecessary parts in the code removed.

EDIT: In response to @LutzL, I attach both the plots of \phi/M_{Pl} and \phi^{'} after integrating forward (solid lines) and backwards (dashed lines), by doing what (s)he said. As you see, there is a sudden deviation from the results of forward integration that I cannot explain.

$|\phi^{'}| vs N$

$\phi/M_{Pl} vs N$

user_3088138 · Answer

I would change the RK4 method to the minimal necessary. It is not necessary to have a v array partially duplicating the contents of the V array, so

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

There are no copy issues here, as v is constructed anew in every step, so that the objects appended to the return array are all separate.

The backward integration should be as simple as the forward integration,

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

and V2[k] should be the same as V1[-k-1] within the bounds of the method error. Large differences are only to be expected in stiff ODE.

Runge-Kutta fourth order method. Integrating backwards

1 Answers