# 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. Commented Jun 11, 2020 at 0:19
• @Wolfgang Bangerth I edited my question to make my problem clear. Commented Jun 11, 2020 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? Commented Jun 11, 2020 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. Commented Jun 11, 2020 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$. Commented Jun 11, 2020 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)$$