Failure of fitting data calculated by Runge-Kutta method to the explicit formulas

Question

I was trying to plot the time-depending function $P(t)$ by Runge-Kutta method. The problem was the R-K fitting data were quite weird from the values calculated by explicit formulas.

The function was solved from a 2nd order ODE of $$\frac{d^2P(t)}{dt^2}=-\frac{t}{\left(t^2+4\right)^{3/2}}$$ with $P(0)=1$ and $\frac{dP(t)}{dt}|_{t=0}=-0.5$

What I did was:

Let $V(t)=\frac{dP(t)}{dt}$, then $A(t)=\frac{d^2P(t)}{dt^2}=\frac{dV(t)}{dt}=-\frac{t}{\left(t^2+4\right)^{3/2}}$,

Let $\mathcal{Y}=\left( \begin{array}{c} \mathcal{Y}_1 \\ \mathcal{Y}_2 \\ \end{array} \right)=\left( \begin{array}{c} V(t) \\ P(t) \\ \end{array} \right)$, then $\frac{d\mathcal{Y}}{dt}=\left( \begin{array}{c} \frac{dV(t)}{dt} \\ \frac{dP(t)}{dt} \\ \end{array} \right)=\left( \begin{array}{c} A(t) \\ V(t) \\ \end{array} \right)=\left( \begin{array}{c} -\frac{t}{\left(t^2+4\right)^{3/2}} \\ \mathcal{Y}_1 \\ \end{array} \right)$.

It's in the form of $\mathcal{Y}'=f(\mathcal{Y},t)$ and $\mathcal{Y}\left(t_0\right)=\left( \begin{array}{c} \mathcal{Y}_1\left(t_0\right) \\ \mathcal{Y}_2\left(t_0\right) \\ \end{array} \right)=\left( \begin{array}{c} -0.5 \\ 1 \\ \end{array} \right)$.

Next, let $t=t_0$ and $w=\mathcal{Y}\left(t_0\right)$, and do iteration with step $h$ $$\left( \begin{array}{c} \text{k1}=h*f(w,t) \\ \text{k2}=h*f\left(\frac{\text{k1}}{2}+w,\frac{h}{2}+t\right) \\ \text{k3}=h*f\left(\frac{\text{k2}}{2}+w,\frac{h}{2}+t\right) \\ \text{k4}=h*f(\text{k3}+w,h+t) \\ w=\frac{1}{6} (\text{k1}+2*\text{k2}+2*\text{k3}+\text{k4})+w \\ t=h+t \\ \end{array} \right)$$

The explicit formulas were easy to worked out as: $$P(t)=-t+\sinh ^{-1}\left(\frac{t}{2}\right)+1$$ However, the R-K fitting seemed to be vibrated as below.

Is there anything wrong in my R-K fitting? Or, is there anything I could do to improve the fitting data?

Thanks.

=======Code=======

Remove["Global`*"];

Y[b_] = ArcSinh[b/2] - b + 1;

ddydx2[y_, t_] = {-(y[[2]]/(4 + (y[[2]])^2)^(3/2)), y[[1]]};

t0 = 0; u0 = {{Y'[t0]}, {Y[t0]}};

h = 0.2; n = Ceiling[40/h];

dx = Table[0, {i, n + 1}];

dy = Table[0, {i, n + 1}, {j, 2}];

w = u0; t = t0;

dx[[1]] = t; dy[[1, ;;]] = Transpose[u0][[1]];

For[i = 1, i <= n, i++, k1 = h * ddydx2[w, t];
  k2 = h * ddydx2[w + k1/2, t + h/2];
  k3 = h * ddydx2[w + k2/2, t + h/2];
  k4 = h * ddydx2[w + k3, t + h];
  w = w + (k1 + 2 * k2 + 2 * k3 + k4)/6; 
  dy[[i + 1, ;;]] = Transpose[w][[1]];
  t = t + h; dx[[i + 1]] = t];

M1 = Transpose[{dx, dy[[;; , 1]]}]; dot1 = ListPlot[M1]; line1 = 
 Plot[Y'[t], {t, t0, t0 + n * h}];

M2 = Transpose[{dx, dy[[;; , 2]]}]; dot2 = ListPlot[M2]; line2 = 
 Plot[Y[t], {t, t0, t0 + n * h}];

Show[dot1, line1]

Show[line2, dot2]

Answer 1

There's a problem with this line

ddydx2[y_, t_] = {-(y[[2]]/(4 + (y[[2]])^2)^(3/2)), y[[1]]};

It should be something like

ddydx2[y_, t_] = {-(t/(4 + (t)^2)^(3/2)), y[[1]]};

This is a version using python

def f(y, t):
    return np.array([y[1], -t / (t**2 + 4)**1.5])

def rk4_step(y, t, h):

    k1 = h * f(y, t)
    k2 = h * f(0.5 * k1 + y, 0.5 * h + t)
    k3 = h * f(0.5 * k2 + y, 0.5 * h + t)
    k4 = h * f(k3 + y, h + t)

    return y + (k1 + 2 * k2 + 2 * k3 + k4) / 6., t + h

def integrate(tmax, h):

    T = [0]
    X = [1]
    V = [-0.5]

    y = np.array([X[0], V[0]])
    t = T[0]
    for k in range(int(tmax / h)):

        y, t = rk4_step(y, t, h)

        X.append(y[0])
        V.append(y[1])
        T.append(t)

    return np.array(T), np.array(X), np.array(V)

Which is obtained by calling

tmax = 10
h = 0.2
t, x, v = integrate(tmax, h)

and then plotting