We slightly modify the example code written to numerically evaluate an integral with the trapezoid rule:
program trapezoidal_integration
implicit none
integer, parameter :: n = 1000 ! Number of trapezoids
real(8) :: a = 0.0d0 ! Lower limit of integration
real(8) :: b = 1.0d0 ! Upper limit of integration
real(8) :: result ! Result of the integration
real(8), dimension(n) :: x_values, f_values
integer :: i
print *, "This program performs numerical integration of f(x) = x^3 from a = ", a, " to b = ", b, " using ", n, " trapezoids."
! Check if n is a valid number of trapezoids
if (n > 0) then
print *, "The number of trapezoids is positive."
else if (n < 0) then
print *, "Error: The number of trapezoids is negative."
stop 1
else
print *, "Error: The number of trapezoids is zero."
stop 1
end if
! Perform numerical integration
result = trapezoidal_rule(a, b, n, x_values, f_values)
! Print the result of the integration
print *, "The integral of f(x) = x^3 from ", a, " to ", b, " is approximately: ", result
! Optionally, print out the x values and their corresponding f(x) values
print *, "x values and f(x) evaluations:"
do i = 1, n-1
print *, "x(", i, ") = ", x_values(i), ", f(x(", i, ")) = ", f_values(i)
end do
contains
! Define the function to integrate: f(x) = x^3
real(8) function f(x)
real(8), intent(in) :: x
f = x**3
end function f
! Trapezoidal rule for numerical integration
real(8) function trapezoidal_rule(a, b, n, x_values, f_values)
real(8), intent(in) :: a, b
integer, intent(in) :: n
real(8), dimension(n), intent(out) :: x_values, f_values
real(8) :: p, sum
integer :: i
p = (b - a) / real(n) ! Width of each trapezoid
sum = 0.5d0 * (f(a) + f(b)) ! End points contribution
! Store the x values and function evaluations in arrays
do i = 1, n-1
x_values(i) = a + i * p
f_values(i) = f(x_values(i)) ! Store the function evaluation
sum = sum + f_values(i)
end do
trapezoidal_rule = sum * p ! Return the computed integral
end function trapezoidal_rule
end program trapezoidal_integration
Two arrays are now introduced in the program:
x_values: This array stores the values where the function is evaluated during the integration.f_values: This array stores the computed values of the function at those values.
In FORTRAN, an array is a collection of elements of the same type stored in contiguous memory locations. Arrays are useful when you need to store multiple values of the same type and access them using an index.
type, dimension(size) :: array_name
where type is the data type of the elements (e.g., integer, real(8), character, etc.), array_name is the name of the array, and size is the number of elements in the array.
[!WARNING]
- Indexing starts at 1: the first element is accessed using
arr(1), notarr(0).- Accessing an index outside the array size (e.g.,
arr(6)when the array has only 5 elements) will typically result in an error.
The function trapezoidal_rule() is modified to take two additional array parameters (x_values and f_values).
real(8) function trapezoidal_rule(func, a, b, n, x_values, f_values)
real(8), intent(in) :: a, b
integer, intent(in) :: n
real(8), dimension(n), intent(out) :: x_values, f_values
...
Inside the do loop, each x value is stored in x_values(i), and the corresponding function evaluation is stored in f_values(i).
do i = 1, n-1
x_values(i) = a + i * p
f_values(i) = f(x_values(i)) ! Store the function evaluation
sum = sum + f_values(i)
end do
After the integration, to print the stored x values and their corresponding function evaluations f(x):
do i = 1, n-1
print *, "x(", i, ") = ", x_values(i), ", f(x(", i, ")) = ", f_values(i)
end do
Array initialisation
Initializing an array means assigning initial values to its elements at the time of declaration. Proper initialization is crucial to avoid unpredictable behavior caused by garbage values.
When you declare an array, you can initialize it in several ways:
- Full Initialization: you provide explicit values for all elements.
integer :: arr(3) = (/1, 2, 3/) ! Array of 3 integers, initialized to {1, 2, 3} - Zero Initialization: You can explicitly initialize all elements to zero (or to any other value):
integer :: arr(4) = 0 ! Array of 4 integers, initialized to {0, 0, 0, 0}
[!WARNING] In FORTRAN, arrays are not automatically initialized when declared. If you declare an array without explicitly initializing it, the array elements will contain garbage values (random values left in memory). This can lead to unpredictable behavior in your program. Remember to always initialise them!
Here an example of arrays declaration and initialisation:
program arrays
implicit none
real, dimension(3) :: a ! Declare a real array 'a' of size 3 (uninitialized)
real :: b(3) ! Declare a real array 'b' of size 3 (uninitialized)
integer, dimension(5) :: c, d ! Declare two integer arrays 'c' and 'd', both of size 5 (uninitialized)
integer :: e(5), f(6) ! Declare two integer arrays 'e' of size 5 and 'f' of size 6 (uninitialized)
real, dimension(3) :: g = (/1.2, 2., 3.7/) ! Declare and initialize a real array 'g' with values 1.2, 2.0, and 3.7
integer, dimension(7) :: h = 1 ! Declare an integer array 'h' of size 7. All elements are initialized with 1
end program
Multidimensional arrays
In Fortran, multidimensional arrays are easily declared by specifying the dimensions in parentheses. These are useful for representing data in matrices, grids, or higher-dimensional structures.
type :: array_name(size1, size2, ...)
An example of a 2D array:
integer :: matrix(3, 4) ! Array of 3 rows and 4 columns
It can be initialised either with a full initialisation:
integer :: matrix(2, 3) = reshape([1, 2, 3, 4, 5, 6], shape(matrix))
Example:
program multidimensional_array
implicit none
integer :: i, j
integer :: matrix(2, 3)
! Initialize the 2D array
matrix = reshape([1, 2, 3, 4, 5, 6], shape(matrix))
! Print the 2D array
print *, "Matrix:"
do i = 1, 2 ! Loop through rows
do j = 1, 3 ! Loop through columns
print *, matrix(i, j),
end do
print *, '' ! Newline after each row
end do
end program multidimensional_array
[!TIP] In Fortran, the
size()function is used to determine the number of elements in an array. It can be applied to arrays of any dimension. The syntax is:n = SIZE(array [, dim])where
arrayis the array whose size you want to determine anddimis optional and representes the dimension along which the size is returned (if omitted,size()returns the total number of elements in the array).! 1D array real :: arr1(5) print *, SIZE(arr1) ! Outputs: 5 ! 2D array real :: arr2(3,4) print *, SIZE(arr2) ! Outputs: 12 (total number of elements) print *, SIZE(arr2, 1) ! Outputs: 3 (size along the first dimension) print *, SIZE(arr2, 2) ! Outputs: 4 (size along the second dimension)
[!WARNING] In Fortran, multidimensional arrays are stored in column-major order, meaning that elements of a column are stored contiguously in memory before moving to the next column. This contrasts with languages like C, which use row-major order, where elements of a row are stored contiguously.
When looping over multidimensional arrays in Fortran, it is more efficient to iterate over the leftmost index (first dimension) in the inner loop and the rightmost index (last dimension) in the outer loop. This is because Fortran stores data in column-major order, so accessing elements in this sequence ensures better memory locality and performance.
program example implicit none real, dimension(128,128,128) :: u ! Velocity field integer :: ix, iy, iz ! Indices used in the do loop do iz = 1, size(u,3) do iy = 1, size(u,2) do ix = 1, size(u,3) ... end do end do end do end programThis approach allows Fortran to access memory contiguously, improving cache performance. If you switch the order of loops (i.e., iterate over the rows in the outer loop and columns in the inner loop), the program will still work but may be less efficient due to non-contiguous memory access.