Two RK4 method in one program

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
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

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.

• What specifically is your problem? It looks like you have a simple system of two differential equations that happen to be coupled. That's totally normal. Books on ODE solvers will tell you how to deal with that. – Wolfgang Bangerth Jun 11 '20 at 0:19
• @Wolfgang Bangerth I edited my question to make my problem clear. – Elham Q Jun 11 '20 at 11:54
• I still don't understand. You have a system of two ODEs. Why do you want to solve them separately, one after the other, when you could solve them together? – Wolfgang Bangerth Jun 11 '20 at 13:54
• @Wolfgang Bangerth Yes you are right but I need a separate code for a. That's why I need to write it in this way. – Elham Q Jun 11 '20 at 14:03
• Then of course you will need to store $a$ at every time where you evaluate it as the right hand side of the ODE for $R$. You will need to put it into an array rather than just a single scalar variable. That's the price you pay: If you want to solve these two ODEs separately, you'll have to store the entire history of $a$ for when you need it for the equation for $R$. – Wolfgang Bangerth Jun 11 '20 at 20:18

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)$$