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}\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):  
        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$

I think, this is also not the correct place to ask. Try to move this to https://scicomp.stackexchange.com — panch, Oct 14 '19 at 13:17
Could you explain deeper the roles of the arrays `V` and `v`? Why is `V` a parameter in the function call? What is the intent behind the line `V[j]=copy.deepcopy(V[j])`? — Lutz Lehmann, Oct 14 '19 at 16:07
@panch : This is not really an algorithm design problem. All the parts of the numerical algorithm are present and correct. What is questionable is the design of the function interface, the assumptions on input and output, which is a programming problem. — Lutz Lehmann, Oct 14 '19 at 16:09
@LutzL: Thanks for commenting. That was part of a more complex code I wrote in which I had to perform several averages and save the results. You can actually ignore it and it won't modify the results here. — , Oct 14 '19 at 16:23
the comments on https://scicomp.stackexchange.com/q/33583/11850 sound sensible, aren't you just trying to do something that isn't valid for that system of equations — Sam Mason, Oct 14 '19 at 16:41
I'm voting to close this question because the scicomp question was answered — Sam Mason, Oct 14 '19 at 17:21

score 0 · Answer 1 · answered Oct 15 '19 at 12:14

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 Answers1