# Runge-Kutta fourth order method. Integrating backwards

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}\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}$$.

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):
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/H)**2)*(1/mpl**2)
return np.array([sigma,-3*(1+HH/3)*sigma-lamb*phi**3/(H**2)])


• Based on my understanding from what you wrote here, I would say your assumption to think that if you integrate an ODE forward in time and reach a final point and then use that final point as initial condition to integrate backward in time and reproduce the same values as in forward integration is simply wrong! Why?! Just look at your ODE... you have $\phi^{'}$ which change its sign with T transformation: $\phi^{'}(-t) = - \phi^{'} (t)$. So, as long as you have that viscous part: $3 (1+\frac{H^{'}}{3H}) \phi^{'}$, nothing will be reversible here... Oct 14, 2019 at 16:01
• You can think of your equation as a Langevin Equation and you should know Langevin Dynamics is not reversible unless there is no viscous or non-conservative force involved which is not the case in your equation. Oct 14, 2019 at 16:03
• @AloneProgrammer: Thanks for commenting. I solved that equation w/o the $\phi^3$ term, and it actually reproduces the values as in forward integration. It seems that term is the responsible, but I don't know why.
– user33042
Oct 14, 2019 at 16:38

Based on my understanding from your question, I assume you are looking for to examine time reversibility of your ODE:

$$\phi^{''}(t) = - 3 (1+\frac{H^{'}(t)}{3 H(t)}) \phi^{'} (t) - \lambda \frac{\phi^{3}(t)}{H^{2}(t)} = 0$$

Where:

$$H(t) = \frac{\lambda \phi^{4}(t)}{4(3M_{Pl}^{2}-0.5\phi^{'}(t)^{2})}$$

and

$$H^{'}(t) = -\frac{1}{2M_{Pl}^{2}} \phi^{'}(t)^{2} H(t)$$

Your initial ODE is in the form of a Langevin equation with a viscous term of $$3 (1+\frac{H^{'}(t)}{3 H(t)}) \phi^{'} (t)$$. Your Langevin equation is not time-reversible due to the fact that the equation will not remain the same under T transformation. Simply because first-order derivative $$\phi^{'}(t)$$ change its sign to negative under T transformation: $$\phi^{'}(-t) = -\phi^{'}(t)$$

• Thanks for commenting, but I also ask why the outcome blow up when going backwards. Thanks to your previous comment, I see why both set of results don't match, but what I don't really understand is why the results diverge in backward integration.
– user33042
Oct 14, 2019 at 16:58
• @J.J It's clear... There is no reason that if your forward integration was stable, the backward one also should be stable... Oct 14, 2019 at 17:00
• @AloneProgrammer One would actually expect it to be unstable if the forward time direction is stable. For example $\dot{x}=A\,x$ is forwards stable if all eigenvalues of $A$ have negative real part. Reversing the time direction would be equivalent to simulating $\dot{x}=-A\,x$ in the forwards direction, which also negates all eigenvalues of $A$ and make it unstable in the backwards direction. Only if there are other equilibria or limit cycles then the solution might stay bounded. Oct 15, 2019 at 13:06