The Project | About Us | Contribute | Donations | License
HOME My courses Completed
◄ Previous
◄

Numerical Methods in Fortran

Next ►
►

Numerical methods are essential for solving mathematical problems that cannot be addressed analytically. Fortran, with its strong numerical computation capabilities, is an excellent language for implementing these methods. In this section, we will cover some fundamental numerical methods and how to implement them in Fortran.

Topics Covered

  1. Introduction to Numerical Methods
  2. Root Finding Methods
  3. Numerical Integration
  4. Solving Ordinary Differential Equations (ODEs)
  5. Practical Exercises

  1. Introduction to Numerical Methods

Numerical methods involve algorithms for approximating solutions to mathematical problems. These methods are crucial in various fields such as engineering, physics, and finance.

Key Concepts

  • Accuracy: How close the numerical solution is to the exact solution.
  • Stability: The behavior of the algorithm when subjected to small changes in input.
  • Efficiency: The computational resources required by the algorithm.

  1. Root Finding Methods

Root finding methods are used to determine the roots of a function, i.e., the values of \( x \) for which \( f(x) = 0 \).

Bisection Method

The bisection method is a simple and robust root-finding method. It works by repeatedly dividing an interval in half and selecting the subinterval that contains the root.

Algorithm

  1. Choose initial interval \([a, b]\) such that \( f(a) \cdot f(b) < 0 \).
  2. Compute the midpoint \( c = \frac{a + b}{2} \).
  3. If \( f(c) = 0 \) or the interval is sufficiently small, return \( c \).
  4. If \( f(a) \cdot f(c) < 0 \), set \( b = c \); otherwise, set \( a = c \).
  5. Repeat steps 2-4 until convergence.

Fortran Implementation

program bisection_method
  implicit none
  real(8) :: a, b, c, tol
  integer :: max_iter, iter

  ! Function prototype
  real(8) :: f(real(8))

  ! Initialize variables
  a = 1.0
  b = 2.0
  tol = 1.0e-6
  max_iter = 100
  iter = 0

  ! Bisection method
  do while (abs(b - a) > tol .and. iter < max_iter)
    c = (a + b) / 2.0
    if (f(c) == 0.0) exit
    if (f(a) * f(c) < 0.0) then
      b = c
    else
      a = c
    end if
    iter = iter + 1
  end do

  ! Output result
  print *, 'Root found at x = ', c
  print *, 'Number of iterations: ', iter

contains
  ! Function definition
  real(8) function f(x)
    real(8), intent(in) :: x
    f = x**3 - x - 2.0
  end function f
end program bisection_method

Explanation

  • Function f: Defines the function whose root we are finding.
  • Variables a and b: Initial interval endpoints.
  • Variable tol: Tolerance for convergence.
  • Variable max_iter: Maximum number of iterations.
  • Loop: Repeatedly halves the interval until the root is found or the maximum number of iterations is reached.

  1. Numerical Integration

Numerical integration approximates the integral of a function over an interval.

Trapezoidal Rule

The trapezoidal rule approximates the integral by dividing the area under the curve into trapezoids.

Formula

\[ \int_a^b f(x) , dx \approx \frac{b - a}{2} \left[ f(a) + f(b) \right] \]

Fortran Implementation

program trapezoidal_rule
  implicit none
  real(8) :: a, b, integral
  integer :: n

  ! Function prototype
  real(8) :: f(real(8))

  ! Initialize variables
  a = 0.0
  b = 1.0
  n = 100

  ! Compute integral
  integral = trapezoidal(a, b, n)

  ! Output result
  print *, 'Integral = ', integral

contains
  ! Function definition
  real(8) function f(x)
    real(8), intent(in) :: x
    f = x**2
  end function f

  ! Trapezoidal rule function
  real(8) function trapezoidal(a, b, n)
    real(8), intent(in) :: a, b
    integer, intent(in) :: n
    real(8) :: h, sum
    integer :: i

    h = (b - a) / n
    sum = 0.5 * (f(a) + f(b))
    do i = 1, n-1
      sum = sum + f(a + i * h)
    end do
    trapezoidal = h * sum
  end function trapezoidal
end program trapezoidal_rule

Explanation

  • Function f: Defines the function to integrate.
  • Variables a and b: Integration limits.
  • Variable n: Number of subintervals.
  • Function trapezoidal: Implements the trapezoidal rule.

  1. Solving Ordinary Differential Equations (ODEs)

ODEs describe the relationship between a function and its derivatives. Numerical methods can approximate solutions to ODEs.

Euler's Method

Euler's method is a simple numerical technique for solving first-order ODEs.

Formula

\[ y_{n+1} = y_n + h f(t_n, y_n) \]

Fortran Implementation

program euler_method
  implicit none
  real(8) :: t, y, h, t_end
  integer :: n, i

  ! Function prototype
  real(8) :: f(real(8), real(8))

  ! Initialize variables
  t = 0.0
  y = 1.0
  h = 0.1
  t_end = 2.0
  n = int((t_end - t) / h)

  ! Euler's method
  do i = 1, n
    y = y + h * f(t, y)
    t = t + h
  end do

  ! Output result
  print *, 'Solution at t = ', t, ' is y = ', y

contains
  ! Function definition
  real(8) function f(t, y)
    real(8), intent(in) :: t, y
    f = -2.0 * t * y
  end function f
end program euler_method

Explanation

  • Function f: Defines the ODE \( \frac{dy}{dt} = f(t, y) \).
  • Variables t and y: Initial conditions.
  • Variable h: Step size.
  • Loop: Iteratively applies Euler's method to approximate the solution.

  1. Practical Exercises

Exercise 1: Implement the Secant Method

The secant method is another root-finding algorithm that uses two initial guesses to approximate the root.

Task

  • Implement the secant method in Fortran.
  • Use the function \( f(x) = x^2 - 4 \).
  • Choose initial guesses \( x_0 = 1.0 \) and \( x_1 = 2.0 \).

Solution

program secant_method
  implicit none
  real(8) :: x0, x1, x2, tol
  integer :: max_iter, iter

  ! Function prototype
  real(8) :: f(real(8))

  ! Initialize variables
  x0 = 1.0
  x1 = 2.0
  tol = 1.0e-6
  max_iter = 100
  iter = 0

  ! Secant method
  do while (abs(x1 - x0) > tol .and. iter < max_iter)
    x2 = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0))
    x0 = x1
    x1 = x2
    iter = iter + 1
  end do

  ! Output result
  print *, 'Root found at x = ', x1
  print *, 'Number of iterations: ', iter

contains
  ! Function definition
  real(8) function f(x)
    real(8), intent(in) :: x
    f = x**2 - 4.0
  end function f
end program secant_method

Exercise 2: Implement Simpson's Rule

Simpson's rule is a more accurate numerical integration method that approximates the integral using parabolic segments.

Task

  • Implement Simpson's rule in Fortran.
  • Use the function \( f(x) = \sin(x) \).
  • Integrate over the interval \([0, \pi]\) with \( n = 100 \).

Solution

program simpsons_rule
  implicit none
  real(8) :: a, b, integral
  integer :: n

  ! Function prototype
  real(8) :: f(real(8))

  ! Initialize variables
  a = 0.0
  b = 3.141592653589793
  n = 100

  ! Compute integral
  integral = simpson(a, b, n)

  ! Output result
  print *, 'Integral = ', integral

contains
  ! Function definition
  real(8) function f(x)
    real(8), intent(in) :: x
    f = sin(x)
  end function f

  ! Simpson's rule function
  real(8) function simpson(a, b, n)
    real(8), intent(in) :: a, b
    integer, intent(in) :: n
    real(8) :: h, sum
    integer :: i

    h = (b - a) / n
    sum = f(a) + f(b)
    do i = 1, n-1
      if (mod(i, 2) == 0) then
        sum = sum + 2.0 * f(a + i * h)
      else
        sum = sum + 4.0 * f(a + i * h)
      end if
    end do
    simpson = h / 3.0 * sum
  end function simpson
end program simpsons_rule

Conclusion

In this section, we explored several fundamental numerical methods and their implementation in Fortran. We covered root-finding methods, numerical integration, and solving ODEs. By practicing these methods, you will gain a deeper understanding of numerical computation and enhance your Fortran programming skills. In the next section, we will delve into more advanced topics and applications of numerical methods in Fortran.

◄ Previous
◄
Next ►
►
The Project | About Us | Contribute | Donations | License
© Copyright 2024. All rights reserved