Speed Comparison:
Fortran vs Julia vs Python
(Runge-Kutta Method)

1. Problem Setting: Newton's Equation of Motion for a One-Dimensional Harmonic Oscillator

To compare the computational speeds of Fortran, Python, and Julia, we consider solving Newton's equation for a one-dimensional harmonic oscillator. The motion equation to be solved is, \[m\frac{d^2}{dt^2}x(t)=-k x(t) \tag{1} \] Thus, Newton's equation of motion is a second-order ordinary differential equation.

To numerically solve Newton's equation, we consider rewriting the second-order ordinary differential equation into the following first-order differential equations. \[\frac{d}{dt}x(t) = v(t) \tag{2}\] \[\frac{d}{dt}v(t) =-\frac{k}{m}x(t) \tag{3} \] The variable \(v(t)\) introduced here is nothing but the velocity. In this analysis, we solve the system of equations (2) and (3) under the initial conditions [\(x(0)=0, v(0)=1\)]. Moreover, we will perform the calculations with both the mass \(m\) and the spring constant \(k\) set to 1.

2. Codes for Fortran, Python, Julia

To compare the computational speeds of Fortran, Python, and Julia, we solve the system of differential equations, equations (2) and (3), using the 4th order Runge-Kutta method. The codes written in each language are shown below. Please refer to the details of the computations in each code. As a rough outline of the code, the time step dt is set to 0.01, and the time evolution calculation is carried out for 10⁸ steps. Additionally, the information of the position \(x(t)\) and velocity \(v(t)\) in the last 1000 steps of the iteration is written to a file.

Fortran Code

(newton.f90)

program main
  implicit none
  real(8) :: x, v
  real(8),allocatable :: xt(:), vt(:)
  real(8) :: mass, k, dt
  integer :: it, nt

  mass = 1d0
  k    = 1d0
  dt   = 1d-2
  nt   = 100000000

  allocate(xt(0:nt), vt(0:nt))

  x = 0d0
  v = 1d0


  do it = 0, nt
    xt(it) = x
    vt(it) = v

    call Runge_Kutta_4th(x,v,dt,mass,k)

  end do

  open(20,file="result_fortran.out")
  do it = nt-1000, nt
    write(20,"(3e26.16e3)")it*dt, xt(it), vt(it)
  end do
  close(20)


  contains 

    subroutine Runge_Kutta_4th(x,v,dt,mass,k)
      implicit none
      real(8),intent(inout) :: x, v
      real(8),intent(in) :: dt, mass, k
      real(8) :: x1,x2,x3,x4,v1,v2,v3,v4

! RK1
      x1 = v
      v1 = force(x, mass, k)

! RK2
      x2 = v+0.5d0*dt*v1
      v2 = force(x+0.5d0*x1*dt, mass, k)

! RK3
      x3 = v+0.5d0*dt*v2
      v3 = force(x+0.5d0*x2*dt, mass, k)

! RK4
      x4 = v+dt*v3
      v4 = force(x+x3*dt, mass, k)

      x = x + (x1+2d0*x2+2d0*x3+x4)*dt/6d0
      v = v + (v1+2d0*v2+2d0*v3+v4)*dt/6d0

    end subroutine Runge_Kutta_4th

    real(8) function force(x,mass,k)
      implicit none
      real(8),intent(in) :: x, mass, k

      force = -x*k/mass

    end function force

end program main

Python code

(newton.py)

import numpy as np

def Runge_Kutta_4th(x, v, dt, mass, k):
    # RK1
    x1 = v
    v1 = force(x, mass, k)

    # RK2
    x2 = v + 0.5 * dt * v1
    v2 = force(x + 0.5 * x1 * dt, mass, k)

    # RK3
    x3 = v + 0.5 * dt * v2
    v3 = force(x + 0.5 * x2 * dt, mass, k)

    # RK4
    x4 = v + dt * v3
    v4 = force(x + x3 * dt, mass, k)

    x = x + (x1 + 2 * x2 + 2 * x3 + x4) * dt / 6
    v = v + (v1 + 2 * v2 + 2 * v3 + v4) * dt / 6

    return x, v

def force(x, mass, k):
    return -x * k / mass

def main():
    mass = 1.0
    k = 1.0
    dt = 1e-2
    nt = 100000000

    xt = np.zeros(nt + 1)
    vt = np.zeros(nt + 1)

    x = 0.0
    v = 1.0

    for it in range(nt + 1):
        xt[it] = x
        vt[it] = v

        x, v = Runge_Kutta_4th(x, v, dt, mass, k)

    with open("result_python.out", "w") as f:
        for it in range(nt - 1000, nt + 1):
            f.write(f"{it * dt:.16e}\t{xt[it]:.16e}\t{vt[it]:.16e}\n")

if __name__ == "__main__":
    main()

Julia code

(newton.jl)

function main()
    mass = 1.0
    k = 1.0
    dt = 1e-2
    nt = 100000000

    xt = zeros(Float64, nt+1)
    vt = zeros(Float64, nt+1)

    x = 0.0
    v = 1.0

    for it = 1:nt+1
        xt[it] = x
        vt[it] = v
        x, v = Runge_Kutta_4th!(x, v, dt, mass, k)
    end

    open("result_julia.out", "w") do file
        for it = nt-999:nt
            println(file, "$(it*dt) $(xt[it]) $(vt[it])")
        end
    end
end

function Runge_Kutta_4th!(x, v, dt, mass, k)
    x1 = v
    v1 = force(x, mass, k)

    x2 = v + 0.5 * dt * v1
    v2 = force(x + 0.5 * x1 * dt, mass, k)

    x3 = v + 0.5 * dt * v2
    v3 = force(x + 0.5 * x2 * dt, mass, k)

    x4 = v + dt * v3
    v4 = force(x + x3 * dt, mass, k)

    x += (x1 + 2 * x2 + 2 * x3 + x4) * dt / 6
    v += (v1 + 2 * v2 + 2 * v3 + v4) * dt / 6

		return x, v
end

function force(x, mass, k)
    return -x * k / mass
end

main()

3. Compilation and Execution of Computations

In this comparison, we use gfortran (GNU Fortran) and ifx (Intel Fortran; formerly ifort) as Fortran compilers. GNU Fortran allows specifying optimization options to improve computational speed at compile time. This time, we will measure the computational speeds for four optimization options: O0, O1, O2, and O3, to also investigate the impact of compilation options on computational speed. For Intel Fortran, in addition to these options, we will add and test the xHOST option, which optimizes according to the executing processor.

The Python code (newton.py) is executed in the terminal with python newton.py, and the Julia code (newton.jl) is similarly executed with julia newton.jl in the terminal. Some may wonder about Julia's compilation time when comparing computational speeds with Fortran, but we have confirmed that this does not affect the results of this comparison, so we will mention this point later.

Executing the prepared Fortran code newton.f90, Python code newton.py, and Julia code newton.jl will output files with the results of solving Newton's equation. The figure below shows the time evolution of the position \(x(t)\) of the harmonic oscillator output by the three codes. It can be seen that all three codes yield identical results.

4. Comparison of Computational Speeds of Fortran, Python, Julia

To compare the computational speeds of Fortran, Python, and Julia, we executed each of the three computation codes five times and measured the computation time for each run. The table below shows the average execution time and standard error obtained from this measurement.

From the table above, let's first compare the results for Fortran compilers and compile options. For Fortran compilers, we used the GNU Fortran compiler (gfortran) and Intel Fortran compiler (ifx; formerly ifort). In the case of the GNU Fortran compiler, specifying the -O0 option to disable optimizations results in execution times about 2.2 times longer than when optimizations are enabled. Therefore, in Fortran, enabling optimizations speeds up computations like solving Newton's equation by about 2.2 times. Also, for GNU Fortran results, there is little difference in execution time when optimization options -O1, -O2, -O3 are used. Thus, for simple computations like this one, higher-level optimizations beyond -O1 do not significantly affect computational speed. Comparing Intel Fortran and GNU Fortran results, execution times for the non-optimized case (-O0) are nearly identical, but computation speeds for Intel Fortran become significantly faster than GNU Fortran when performing optimizations of -O1 or higher. Especially, performing -O3 optimization along with -xHOST (optimization according to the processor) makes Intel Fortran's computation speed about 1.8 times faster than GNU Fortran.

Comparing the execution time of Fortran (ifx -O3 -xHOST) with Python, Fortran is about 81.9 times faster than Python. Moreover, comparing Fortran (ifx -O3 -xHOST) with Julia, Fortran is about 3.3 times faster than Julia. Furthermore, comparing the execution time of Julia with Python, Julia is about 24.6 times faster than Python.

[Back to Research and Verification Home]