0
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – Wolfgang Bangerth Jun 11 at 0:19
  • $\begingroup$ @Wolfgang Bangerth I edited my question to make my problem clear. $\endgroup$ – Elham Q Jun 11 at 11:54
  • $\begingroup$ 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? $\endgroup$ – Wolfgang Bangerth Jun 11 at 13:54
  • $\begingroup$ @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. $\endgroup$ – Elham Q Jun 11 at 14:03
  • 2
    $\begingroup$ 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$. $\endgroup$ – Wolfgang Bangerth Jun 11 at 20:18
1
$\begingroup$

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

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.