Magnetostatics in Mathematica and Fortran

Aleksandar Donev - Dr. Phillip Duxbury

November, 2000

Help with Fortran Solution

Trapezoid-Rule for Numerical Integration

It was purely fortunate (due to the relatively simple geometry) that the integrals above could be evaluated analytically. In most cases this has do be done numerically. Fortran libraries by far surpass Mathematica in speed when it comes to such numerical manipulations, but they are not nearly as easy to use.

Before we attack the integration problem in eq. , we need to generalize it to a standard problem that can be applied to a variety of problems and that is well studied in numerical analysis,

$\begin{displaymath}I=\int_{a}^{b}f(x)dx \end{displaymath}$

(1)

where f(x) is more-or-less a well behaved function of x in [a,b]. There are many so-called quadrature formulae for numerically calculating I, such as Newton-Cotes, Gauss and Romberg quadrature. A nice desktop reference for these and many other numerical algorithms is ``Numerical Recipies in FORTRAN 77 and Fortran 90'', found at, for example, http://lib-www.lanl.gov/numerical/.

The simplest, yet robust, quadrature is the trapezoid-rule integration. It starts by splitting the interval [a,b] into N intervals of equal length $h=\frac{\left\vert a-b\right\vert }{N}$ with positions x_i=a+ih , 0<i<N. Then it is true that,

$\begin{displaymath}\int_{a}^{b}f(x)dx\approx h\sum_{i=1}^{N-1}f(x_{i})+\frac{h}{2}[f(a)+f(b)% ]+O(h^{2}) \end{displaymath}$

(2)

where O(h²) denotes error of order h². It is easy to see that it is not at all difficult to code in Fortran. Write a function, say Trapezoid, that integrates a given function f(x) on [a,b]. We've had examples of passing procedures as arguments to procedures in the graphics routines, and although this is not in the manual, it's not hard to understand. All one needs to do is provide an interface for the procedure:

   FUNCTION Trapezoid(f,a,b,N) RESULT(I)
      INTERFACE
         FUNCTION f(x) RESULT(f_x) ! Function to be integrated
           REAL, INTENT(IN) :: x
           REAL :: f_x
         END FUNCTION f                   
      END INTERFACE        
      REAL, INTENT(IN) :: a,b ! The bounds
      INTEGER, INTENT(IN) :: N
      ...
      ...Calculate I here using trapezoid rule using f and a DO loop...
      ...
   END FUNCTION Trapezoid

Test this procedure on some simple example you already know the analytical answer for. Next we describe a ready-made quadrature library. You don't have to use it, but it isn't much different from using Trapezoid.

Advanced: Module `Numerical_Integration`

As usual, we made a simple wrapper routine for numerical integration, AdaptiveIntegral, which in turn calls a ready-made adaptive quadrature Fortran library taken from Adam Miller's Fortran 90 website. The function is in the module Numerical_Integration, which is in our class directory and you must USE in your program. This function accepts a user function f to be integrated in a given finite range [a,b] and returns the integral to within the desired relative or absolute precisions relative_error and absolute_error (these are both optional, so you don't need to specify them). An estimate of the error in the result is returned in the argument error_estimate, and the function also returns an error signal if the integration failed (this is placed in the module variable error_flag and you need not worry about it)

The interface to this function and the module variables are given here:

module Numerical_Integration
  use Precision, only: wp ! wp=dp in the module Precision
  use Adapt_Quad ! Read-made library
  private
 
  public :: AdaptiveIntegral ! The main routine
  integer, save, public :: error_flag=0 ! An error flag (0 is success)
 
  contains
  function AdaptiveIntegral(f,a,b,error_estimate, &
                            relative_error,absolute_error) &
                            result(integral)
 
    ! f is the function we want to integrate:
    interface
      function f(x) result(f_x)
        use Precision, only: wp
        real(kind=wp), intent(in) :: x
        real(kind=wp) :: f_x
      end function f
    end interface  
    real(kind=wp), intent(in) :: a,b ! The interval of integration
    real(kind=wp), intent(out) :: error_estimate ! Estimate of error
    ! The desired errors (they are optional): 
    real(kind=wp), intent(in), optional :: relative_error,absolute_error  
    real(kind=wp) :: integral ! The result of the integration
    ...
  end function AdaptiveIntegral   
 
end module Numerical_Integration

It is important to point that in the module Precision it is defined that wp=dp, so you need to use double precision in your main program as well (this library only came in double precision). There will be no examples given of how to call this routine. You should be able to read the above specifications and call this function without any problems.

The Module `B_theta`

Just like the plotting routines in the FunGraphics module, the above integration routine needs to be supplied with a user function to integrate. In our case, this function will be related to $\vec{B}_{\theta }$ , since this is the integrand we need to integrate. Use the expression that Mathematica gave you for $\vec{B}_{\theta }$ to write two functions, say dBx and dBy, that have the same interface as the function f in the routine AdaptiveIntegral (see above) and return the x and z components of $\vec{B}_{\theta }$ respectively.

These will accept $\theta$ as an input argument (why?), but will need to know the values of x and z. By placing the routines in a module, say B_theta, and making x and z global public variables in this module, the main program can set these values (say the user enters them) before calling the routine AdaptiveIntegral, and all should work well. Re-read the description of this methodology given in the last worksheet so you understand how it works.

Now compile and run your program and evaluate the magnetic field at the same point as you did in your Mathematica worksheet. Compare the two results.

Plotting the Magnetic Field in Fortran

Once you make the above work, if you have time and will, you can plot the magnetic field from within Fortran using the FunContourPlot2D routine from the last worksheet. If you want some Fortran challenge, do this and we will help you. One thing to carry in mind is that the plotting routines accept only single precision arguments, while the integration routines above accept only double precision variables. So you need to be more careful this time than the usual REAL(KIND=wp) in all declarations.

You will get extra credit applied toward your final exam if you do this!

Compiling the Program

As usual, the ready-made files are in our class directory. Before you compile your code, execute:

#   source /classes/phy201/Integration.init

and then append a $ADQDvf90 (stands for adaptive quadrature with VAST f90) to every compilation line. If you use the graphics modules as well source the file SimpleGraphics.init as well and append a $DISLINvf90 as well.

About this document ...

Magnetostatics in Mathematica and Fortran

This document was generated using the LaTeX2HTML translator Version 98.1p1 release (March 2nd, 1998)

The command line arguments were:
latex2html -split 0 worksheet8_fhelp.tex.

The translation was initiated by Phil Duxbury on 2000-11-06

Phil Duxbury
2000-11-06