Two RK4 method in one program

Question

I want to solve this integral using RK4 by coding in Fortran:

$$R=∫1/a(t) dt → dR/dt=1/a(t) =f(t)$$

Initial point: t=0 (or a=0.001) and R=0

And I have to get a(t) by solving another differential equation: $$da/dt=1/a+1/a^2 =g(a)$$

Initial point: t=0 and a=0.001

I wrote this code to get a(t):

    PROGRAM RK4
      implicit none
      real h,t
      integer n
      read*,h,n
      call Scale_Factor(h,n,t,a)
    END PROGRAM

    !---------------------------------------------

    SUBROUTINE Scale_Factor(h,n,t,a)
      implicit none
      real t,a,k1,k2,k3,k4,h,g
      integer i,n

      t=0
      a=0.001


    Do i=1,n

       k1=h*g(a)

       k2=h*g(a+k1/2.0)

       k3=h*g(a+k2/2.0)

       k4=h*g(a+k3)

       t=t+h

       a=a+(k1+2*k2+2*k3+k4)*(1/6.0)

       write(*,*)t,a

    END DO
    END SUBROUTINE

    !-------------------------
    FUNCTION g(a)
      implicit none
      real a,g
      g=sqrt((1.0/a)+(1.0/a**2)) 
    END FUNCTION

And I have another similar program for solving the first integral. But I need to use a(t) that this program produces to solve the integral and I do not know how to combine them in a single program.

What I wrote to combine them is this:

    Program RK4
    implicit none

    real k1,k2,k3,k4,h,t,R
    integer i,n
    real a
    read*,n,h 


    t=0
    R=0

    Do i=1,n

      k1=h*(1/a(t))

      k2=h*(1/a(t+h/2.0))

      k3=h*(1/a(t+h/2.0))

      k4=h*(1/a(t+h))

      t=t+h

      R=R+(k1+2*k2+2*k3+k4)*(1/6.0)

      write(*,*)t,R

    End Do

    end program

    !-----------------------------------------

    SUBROUTINE Scale_Factor(h,n,t,a)
      implicit none
      real t,a,k1,k2,k3,k4,h,g
      integer i,n

      t=0
      a=0.001


    Do i=1,n

       k1=h*g(a)

       k2=h*g(a+k1/2.0)

       k3=h*g(a+k2/2.0)

       k4=h*g(a+k3)

       t=t+h

       a=a+(k1+2*k2+2*k3+k4)*(1/6.0)

       write(*,*)t,a

    END DO
    END SUBROUTINE

    !-------------------------
    FUNCTION g(a)
      implicit none
      real a,g
      g=sqrt((1.0/a)+(1.0/a**2)) 
    END FUNCTION

But I know it is not correct.

Answer 1

First one should verify the numerical solution for $a(t)$ against the analytic solution.

$ d_t a = 1/a + 1/a^2 = \frac{a+1}{a^2} $

Use $b=a+1$, then the ODE becomes

$ \frac{(b-1)^2}{b} db = dt, $

from which we find the solution

$ b^2/2 - 2 b - \ln(b) = t $

This relation (defined up to a constant to satisfy the initial conditions) is not invertible but it is sufficient for verifying the numerical solution for $a(t) = b(t)-1$.

Once $a(t)$ is available, it is trivial to solve numerically the other ODE,

$ d_t R = 1/a $

To verify $R(t)$ one can find it analytically as well,

$ \frac{dR}{dt} = \frac{dR}{da} \frac{da}{dt} = \frac{1}{a} $

Using $b=a+1$ this can be rewritten as

$ dR = \frac{b-1}{b} db, $

which leads to the solution (up to an additive constant to satisfy initial conditions)

$ R = b - \ln(b) $