Fortran scientific simulations
last modified November 14, 2026
A simulation replaces a physical process by a rule that is applied again and again: a cellular automaton updates a grid of cells, a differential equation advances a set of numbers by one time step, a random walk adds a random displacement to a position. The code is usually short, and the difficulty is not in the loops but in the choices around them: how large a time step keeps the method stable, which quantity should stay constant, and how the result should be stored.
This article writes seven small simulations in Fortran: the Game of Life, the one-dimensional heat equation, a planet in a circular orbit, a two-dimensional random walk, the same heat equation saved as a PGM image, a pendulum integrated by two methods of different quality, and the predator-prey equations. Each one is checked against something that must not change, an energy, an angular momentum, an invariant or a known mean.
The Game of Life
The Game of Life is a cellular automaton: a grid of cells, each alive or dead, updated in rounds. A live cell survives with two or three live neighbours and dies otherwise; a dead cell with exactly three live neighbours becomes alive; every cell is updated from the previous round at the same time, so the program needs two grids, the current one and the next one. The pattern below is a glider, five live cells that walk across the grid.
program k1
implicit none
integer, parameter :: n = 8, rounds = 6
logical :: grid(n, n), next(n, n)
integer :: i, j, k, neighbours
grid = .false.
grid(2, 3) = .true.
grid(3, 4) = .true.
grid(4, 2) = .true.
grid(4, 3) = .true.
grid(4, 4) = .true.
do k = 0, rounds
print '(a, i0, a, i0, a)', 'generation ', k, ', population ', &
count(grid), ':'
do i = 1, n
print '(8a1)', (merge('#', '.', grid(i, j)), j = 1, n)
end do
if (k == rounds) exit
next = .false.
do i = 2, n - 1
do j = 2, n - 1
neighbours = count(grid(i - 1:i + 1, j - 1:j + 1)) - &
merge(1, 0, grid(i, j))
next(i, j) = neighbours == 3 .or. &
(neighbours == 2 .and. grid(i, j))
end do
end do
grid = next
end do
end program k1
The neighbour count uses an array section: grid(i-1:i+1, j-1:j+1)
is the three-by-three block around a cell, count returns how many of
its nine cells are alive, and subtracting the cell itself gives the eight
neighbours. The boundaries of the array act as a dead border because the loops
stop one cell short of the edge, which is enough for a glider in the middle. The
five cells keep their number and move one cell diagonally every four rounds,
which the six printed generations show step by step.
$ gfortran k1.f90 -o k1 $ ./k1 generation 0, population 5: ........ ..#..... ...#.... .###.... ........ ........ ........ ........ generation 1, population 5: ........ ........ .#.#.... ..##.... ..#..... ........ ........ ........ generation 2, population 5: ........ ........ ...#.... .#.#.... ..##.... ........ ........ ........ generation 3, population 5: ........ ........ ..#..... ...##... ..##.... ........ ........ ........ generation 4, population 5: ........ ........ ...#.... ....#... ..###... ........ ........ ........ generation 5, population 5: ........ ........ ........ ..#.#... ...##... ...#.... ........ ........ generation 6, population 5: ........ ........ ........ ....#... ..#.#... ...##... ........ ........
The heat equation in one dimension
The temperature along a thin rod obeys the heat equation, which says that a
point warms in proportion to how much hotter its neighbours are. Discretised on
a grid of points it becomes an update rule: the new value at a point is the old
value plus a fraction of the difference between the average of its two
neighbours and itself. The fraction, written r below, combines the
conductivity, the time step and the spacing, and the scheme is stable only while
it stays at or below one half.
program k2
implicit none
integer, parameter :: n = 11, steps = 100
real, parameter :: r = 0.4
real :: u(n), following(n)
integer :: i, k
u = 0.0
u(1) = 100.0
print '(a)', 'step temperature at eleven points of the rod'
call profile(0, u)
do k = 1, steps
following(1) = 100.0
following(n) = 0.0
do i = 2, n - 1
following(i) = u(i) + r * (u(i - 1) - 2.0 * u(i) + u(i + 1))
end do
u = following
if (mod(k, 20) == 0) call profile(k, u)
end do
contains
subroutine profile(step, values)
integer, intent(in) :: step
real, intent(in) :: values(:)
print '(i4, 11f6.1)', step, values
end subroutine profile
end program k2
The rod starts at zero everywhere except the left end, which is held at 100
degrees, and the right end, which is held at zero. Heat spreads inwards and the
profile becomes smoother with every row, until it approaches the straight line
from 100 to 0 that carries a constant flux. The value 0.4 for r
keeps the scheme stable; a value above 0.5 would make every step amplify the
rounding errors and the numbers would oscillate and grow without bound.
$ gfortran k2.f90 -o k2 $ ./k2 step temperature at eleven points of the rod 0 100.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 20 100.0 80.6 62.2 46.0 32.3 21.6 13.6 8.1 4.4 1.9 0.0 40 100.0 86.0 72.5 59.6 47.8 37.2 27.9 19.7 12.5 6.1 0.0 60 100.0 88.2 76.6 65.3 54.5 44.3 34.5 25.4 16.6 8.2 0.0 80 100.0 89.2 78.5 67.9 57.5 47.4 37.5 27.9 18.5 9.2 0.0 100 100.0 89.6 79.3 69.1 58.9 48.8 38.9 29.1 19.3 9.6 0.0
An orbit that keeps its energy
A planet around a star obeys Newton's law with an inverse-square force, and a circular orbit is a good test of an integrator because its exact answer is known: the radius stays one, the period is two pi, and both the energy and the angular momentum are constants of the motion. The program uses velocity Verlet, which advances the position by a half step of velocity, applies the acceleration, and finishes the velocity step; it is one of the family of symplectic methods, built so that the energy oscillates around the true value instead of drifting away from it.
program k3
use iso_fortran_env, only: real64
implicit none
integer, parameter :: steps = 6280
real(real64), parameter :: dt = 0.001_real64
real(real64) :: r(2), v(2), a(2), energy, initial, angular
integer :: k
r = [1.0_real64, 0.0_real64]
v = [0.0_real64, 1.0_real64]
initial = energy_of(r, v)
print '(a)', ' t x y energy drift angular'
do k = 1, steps
call accelerate(r, a)
v = v + 0.5_real64 * dt * a
r = r + dt * v
call accelerate(r, a)
v = v + 0.5_real64 * dt * a
if (mod(k, 628) == 0) then
energy = energy_of(r, v)
angular = r(1) * v(2) - r(2) * v(1)
print '(f6.3, 2f9.5, es12.4, es10.2, f9.5)', real(k, real64) * &
dt, r, energy, abs(energy - initial), angular
end if
end do
print '(a, f0.6)', 'a full circle takes 2*pi = ', 2.0_real64 * &
3.141592653589793_real64
print '(a, f0.6)', 'the radius after one period: ', norm2(r)
contains
subroutine accelerate(position, acceleration)
real(real64), intent(in) :: position(2)
real(real64), intent(out) :: acceleration(2)
acceleration = -position / norm2(position)**3
end subroutine accelerate
real(real64) function energy_of(position, velocity) result(e)
real(real64), intent(in) :: position(2), velocity(2)
e = 0.5_real64 * dot_product(velocity, velocity) - &
1.0_real64 / norm2(position)
end function energy_of
end program k3
The planet passes through the four cardinal points on its way around and returns to a radius of 1.000000 after 6.28 time units, which is two pi to the printed digits. The energy is -0.5 throughout, as it must be for a circular orbit of radius one with these units, and its drift never exceeds 1.3e-13. The angular momentum is exactly 1.00000 in every row, because the force always points at the centre and cannot change it. A method that did not preserve these quantities would spiral inwards or outwards within a few orbits.
$ gfortran k3.f90 -o k3
$ ./k3
t x y energy drift angular
0.628 0.80920 0.58753 -5.0000E-01 1.05E-14 1.00000
1.256 0.30962 0.95086 -5.0000E-01 4.32E-14 1.00000
1.884 -0.30811 0.95135 -5.0000E-01 8.17E-14 1.00000
2.512 -0.80827 0.58882 -5.0000E-01 1.14E-13 1.00000
3.140 -1.00000 0.00159 -5.0000E-01 1.27E-13 1.00000
3.768 -0.81014 -0.58624 -5.0000E-01 1.16E-13 1.00000
4.396 -0.31114 -0.95036 -5.0000E-01 8.45E-14 1.00000
5.024 0.30659 -0.95184 -5.0000E-01 4.65E-14 1.00000
5.652 0.80733 -0.59010 -5.0000E-01 1.44E-14 1.00000
6.280 0.99999 -0.00319 -5.0000E-01 3.05E-15 1.00000
a full circle takes 2*pi = 6.283185
the radius after one period: 1.000000
A random walk and its diffusion
A random walk takes a fixed number of steps, each in one of four directions with
equal probability, and asks how far from the start the walker ends up. The
answer is not a single number but a distribution, and its mean square grows
linearly with the number of steps: n steps of length one give a mean
square distance of n, which is the discrete version of diffusion.
Averaging four hundred walks at four lengths makes the law visible.
program k4
use iso_fortran_env, only: real64
implicit none
integer, parameter :: trials = 400
integer, parameter :: lengths(4) = [100, 400, 1600, 6400]
integer :: i, t, trial, steps, direction
real(real64) :: x, y, total, u
call seed_with(4242)
print '(a)', ' steps mean square distance plain steps'
do t = 1, size(lengths)
steps = lengths(t)
total = 0.0_real64
do trial = 1, trials
x = 0.0_real64
y = 0.0_real64
do i = 1, steps
call random_number(u)
direction = int(u * 4.0_real64)
select case (direction)
case (0)
x = x + 1.0_real64
case (1)
x = x - 1.0_real64
case (2)
y = y + 1.0_real64
case (3)
y = y - 1.0_real64
end select
end do
total = total + x * x + y * y
end do
print '(i7, f22.1, f14.0)', steps, total / trials, real(steps, real64)
end do
contains
subroutine seed_with(base)
integer, intent(in) :: base
integer, allocatable :: seed(:)
integer :: i, n
call random_seed(size=n)
allocate (seed(n))
do i = 1, n
seed(i) = base + 37 * i
end do
call random_seed(put=seed)
end subroutine seed_with
end program k4
The measured mean square distances are 101, 382, 1577 and 6584 against the predicted 100, 400, 1600 and 6400. Each estimate is an average over four hundred walks, so it carries a sampling error of a few percent, and the errors do not shrink with the length of the walk, only with the number of walks. The seed is fixed at the start, so the table is the same on every run, which makes it possible to compare a change in the program with the same random numbers.
$ gfortran k4.f90 -o k4
$ ./k4
steps mean square distance plain steps
100 101.0 100.
400 381.6 400.
1600 1577.1 1600.
6400 6583.5 6400.
Heat as a grayscale image
A two-dimensional simulation produces a picture, and the simplest format for one
is PGM, the portable graymap: a short text header that gives the magic
number P5, the width, the height and the maximum value, followed by
one byte per pixel in row order. The program below diffuses a hot square on a
64 by 64 grid and writes the temperature as brightness, using stream access so
that the bytes land in the file exactly as they are stored in the string.
program k5
implicit none
integer, parameter :: n = 64, steps = 80
real, parameter :: r = 0.2
real :: u(n, n), following(n, n)
character(len=n * n) :: pixels
character(len=64) :: header
character(len=:), allocatable :: image
integer :: i, j, k, unit, bytes, low, high, value, position
u = 0.0
do i = 26, 38
do j = 26, 38
u(i, j) = 255.0
end do
end do
do k = 1, steps
following = u
do i = 2, n - 1
do j = 2, n - 1
following(i, j) = u(i, j) + r * (u(i - 1, j) + u(i + 1, j) + &
u(i, j - 1) + u(i, j + 1) - 4.0 * u(i, j))
end do
end do
u = following
end do
position = 0
do i = 1, n
do j = 1, n
position = position + 1
value = min(255, max(0, nint(u(i, j))))
pixels(position:position) = achar(value)
end do
end do
write (header, '(a, 3(1x, i0))') 'P5', n, n, 255
open (newunit=unit, file='heat.pgm', access='stream', &
form='unformatted', status='replace')
write (unit) trim(header)//char(10)
write (unit) pixels
close (unit)
inquire (file='heat.pgm', size=bytes)
allocate (character(len=bytes) :: image)
open (newunit=unit, file='heat.pgm', access='stream', &
form='unformatted', status='old')
read (unit) image
close (unit)
low = 255
high = 0
do i = len_trim(header) + 2, bytes
value = iachar(image(i:i))
low = min(low, value)
high = max(high, value)
end do
print '(a, i0)', 'file size in bytes: ', bytes
print '(a, a)', 'the header reads: ', image(1:len_trim(header))
print '(a, i0)', 'pixel bytes: ', bytes - len_trim(header) - 1
print '(a, i0, a, i0)', 'the darkest pixel is ', low, &
', the brightest ', high
print '(a)', 'a coarse view of the temperature:'
do i = 1, n, 4
do j = 1, n, 4
value = nint(u(i, j))
if (value * 4 > high * 3) then
write (*, '(a)', advance='no') '#'
else if (value * 2 > high) then
write (*, '(a)', advance='no') '+'
else if (value * 10 > high) then
write (*, '(a)', advance='no') '.'
else
write (*, '(a)', advance='no') ' '
end if
end do
print *
end do
end program k5
The file is 4109 bytes: thirteen for the header and its newline, and 4096 for the pixels. The program reads the file back, checks its size, takes the header apart and scans the bytes for the darkest and the brightest pixel, which are 0 and 143 after eighty rounds. The coarse text view at the end samples every fourth cell in both directions, so each character stands for a four-by-four block. Any image viewer opens the file, and the same technique writes the frames of a simulation as a numbered sequence.
$ gfortran k5.f90 -o k5
$ ./k5
file size in bytes: 4109
the header reads: P5 64 64 255
pixel bytes: 4096
the darkest pixel is 0, the brightest 143
a coarse view of the temperature:
.....
..++...
.+##+..
.+###..
..+#+..
......
....
A pendulum, integrated twice
The pendulum equation theta'' = -(g/L) sin(theta) has no closed
solution in elementary functions, so it is a natural test case. The program
integrates it twice from the same starting angle, once with the classical
fourth-order Runge-Kutta method and once with the simplest explicit Euler step,
and compares the total energy of each with the energy it started from.
Runge-Kutta evaluates the slope four times per step; Euler evaluates it once and
is first-order accurate.
program k6
use iso_fortran_env, only: real64
implicit none
real(real64), parameter :: g = 9.81_real64, length = 1.0_real64, &
dt = 0.01_real64
real(real64) :: theta, omega, theta2, omega2, e0
integer :: k
theta = 1.0_real64
omega = 0.0_real64
theta2 = theta
omega2 = omega
e0 = energy(theta, omega)
print '(a)', ' t rk4 drift euler drift rk4 angle euler angle'
do k = 1, 1200
call rk4_step(theta, omega, dt)
call euler_step(theta2, omega2, dt)
if (mod(k, 300) == 0) then
print '(f6.2, 2es14.2, 2f12.6)', real(k, real64) * dt, &
energy(theta, omega) - e0, energy(theta2, omega2) - e0, &
theta, theta2
end if
end do
contains
real(real64) function acceleration(angle) result(a)
real(real64), intent(in) :: angle
a = -(g / length) * sin(angle)
end function acceleration
real(real64) function energy(angle, rate) result(e)
real(real64), intent(in) :: angle, rate
e = 0.5_real64 * length * rate * rate + g * (1.0_real64 - cos(angle))
end function energy
subroutine rk4_step(angle, rate, h)
real(real64), intent(inout) :: angle, rate
real(real64), intent(in) :: h
real(real64) :: a1, a2, a3, a4, v1, v2, v3, v4
v1 = rate
a1 = acceleration(angle)
v2 = rate + 0.5_real64 * h * a1
a2 = acceleration(angle + 0.5_real64 * h * v1)
v3 = rate + 0.5_real64 * h * a2
a3 = acceleration(angle + 0.5_real64 * h * v2)
v4 = rate + h * a3
a4 = acceleration(angle + h * v3)
angle = angle + h * (v1 + 2.0_real64 * v2 + 2.0_real64 * v3 + v4) / &
6.0_real64
rate = rate + h * (a1 + 2.0_real64 * a2 + 2.0_real64 * a3 + a4) / &
6.0_real64
end subroutine rk4_step
subroutine euler_step(angle, rate, h)
real(real64), intent(inout) :: angle, rate
real(real64), intent(in) :: h
rate = rate + h * acceleration(angle)
angle = angle + h * rate
end subroutine euler_step
end program k6
The energy of the Runge-Kutta solution wanders by a few times 1e-8 over twelve seconds, while the Euler solution is off by up to 6e-2, more than a percent of the total, and its angle is wrong in the third decimal. Neither drift is a rounding error: both are the truncation error of the method, and the fourth-order method reduces it by a factor of about a million at this step size. The pendulum swings between one radian and minus one, as it should, since it starts at rest from that angle.
$ gfortran k6.f90 -o k6 $ ./k6 t rk4 drift euler drift rk4 angle euler angle 3.00 -2.02E-08 5.99E-02 -0.823834 -0.832202 6.00 -2.58E-08 4.82E-02 0.344653 0.359047 9.00 -3.58E-08 -3.59E-02 0.270586 0.255599 12.00 -5.82E-08 -6.27E-02 -0.778413 -0.768739
Predator and prey
The Lotka-Volterra equations describe two populations, prey that grow without limit and predators that eat them: the prey grow at a rate proportional to their number, lose some to the predators, and the predators grow in proportion to the prey they catch and die at a constant rate. The system has a conserved quantity, a combination of the two populations that must stay constant along an exact trajectory, which makes it a strict test of the integration.
program k7
use iso_fortran_env, only: real64
implicit none
real(real64), parameter :: alpha = 1.1_real64, beta = 0.4_real64, &
gamma = 0.1_real64, delta = 0.4_real64, dt = 0.01_real64
real(real64) :: prey, predator, start, drift_max, low, high
integer :: k
prey = 10.0_real64
predator = 5.0_real64
start = invariant(prey, predator)
drift_max = 0.0_real64
low = prey
high = prey
print '(a)', ' t prey predator drift of the invariant'
do k = 1, 4000
call rk4_step(prey, predator, dt)
drift_max = max(drift_max, abs(invariant(prey, predator) - start))
low = min(low, prey)
high = max(high, prey)
if (mod(k, 400) == 0) then
print '(f6.1, 2f13.5, es16.2)', real(k, real64) * dt, prey, &
predator, invariant(prey, predator) - start
end if
end do
print '(a, f0.3, a, f0.3)', 'prey runs from ', low, ' to ', high
print '(a, es10.2)', 'the largest drift of the invariant: ', drift_max
print '(a, 2f9.3)', 'the equilibrium point: ', delta / gamma, &
alpha / beta
contains
pure real(real64) function growth(x, y) result(r)
real(real64), intent(in) :: x, y
r = alpha * x - beta * x * y
end function growth
pure real(real64) function loss(x, y) result(r)
real(real64), intent(in) :: x, y
r = gamma * x * y - delta * y
end function loss
pure real(real64) function invariant(x, y) result(v)
real(real64), intent(in) :: x, y
v = gamma * x - delta * log(x) + beta * y - alpha * log(y)
end function invariant
subroutine rk4_step(x, y, h)
real(real64), intent(inout) :: x, y
real(real64), intent(in) :: h
real(real64) :: k1x, k1y, k2x, k2y, k3x, k3y, k4x, k4y
k1x = growth(x, y)
k1y = loss(x, y)
k2x = growth(x + 0.5_real64 * h * k1x, y + 0.5_real64 * h * k1y)
k2y = loss(x + 0.5_real64 * h * k1x, y + 0.5_real64 * h * k1y)
k3x = growth(x + 0.5_real64 * h * k2x, y + 0.5_real64 * h * k2y)
k3y = loss(x + 0.5_real64 * h * k2x, y + 0.5_real64 * h * k2y)
k4x = growth(x + h * k3x, y + h * k3y)
k4y = loss(x + h * k3x, y + h * k3y)
x = x + h * (k1x + 2.0_real64 * k2x + 2.0_real64 * k3x + k4x) / &
6.0_real64
y = y + h * (k1y + 2.0_real64 * k2y + 2.0_real64 * k3y + k4y) / &
6.0_real64
end subroutine rk4_step
end program k7
The populations cycle: the prey run up to 13.7 while the predators are few, then the predators multiply and the prey fall to 0.5, and the cycle begins again. The invariant stays constant to 3e-10 over forty time units, which is the accuracy of the fourth-order method at this step. The equilibrium point, where births and deaths balance, is at four prey and 2.75 predators, computed from the four coefficients at the end. A run with a larger step would show the cycle slowly drifting away from its starting loop.
$ gfortran k7.f90 -o k7 $ ./k7 t prey predator drift of the invariant 4.0 0.51944 2.44324 -2.94E-11 8.0 4.48460 0.91917 -2.58E-11 12.0 2.19315 5.88055 -8.99E-11 16.0 0.74741 1.61980 -1.18E-10 20.0 10.21411 1.33381 -1.40E-10 24.0 0.71348 4.20571 -2.19E-10 28.0 1.47109 1.13294 -2.07E-10 32.0 12.75553 3.79926 -2.34E-10 36.0 0.50860 2.77656 -3.05E-10 40.0 3.42768 0.92145 -2.97E-10 prey runs from .509 to 13.676 the largest drift of the invariant: 3.11E-10 the equilibrium point: 4.000 2.750
Next steps
The arrays article explains the two-dimensional
grids and the array sections that take the neighbour counts, and the
binary files article covers the stream
access used to write the image. The
random numbers article describes the
seeding that makes the walk reproducible, and the
numerical integration article
compares the accuracy of the same methods on a scalar problem. The
loops article covers the counted and conditional
loops of every simulation, and the
formatted output article explains the
es and f descriptors that print the tables.
Source
Fortran 2018 standard (J3/18-007r1)
In this article we have run the Game of Life on a grid with two generations in
memory, solved the one-dimensional heat equation with an explicit scheme and
watched the profile relax to a straight line, integrated a circular orbit with
velocity Verlet and checked the energy and the angular momentum, measured the
diffusion of a random walk against the n law, written a two-
dimensional temperature field as a PGM image and read it back, compared
Runge-Kutta with Euler on a pendulum, and followed the predator-prey cycle while
checking its invariant. The previous article,
Fortran linear algebra, supplies the
matrix operations that a larger simulation would need.
Author
List all Fortran tutorials.