Fortran Example – Mathieu functions and characteristics¶

Here we consider a set of functions that can not (generally) be written in closed form. The Mathieu functions, \(y(x)\), are defined as the solutions of the differential equation:

(1)¶\[\frac{d^2y}{dx^2} + \left(a-2q\cos(2x)\right)y(x) = 0\]

where \(a,q\in\mathbb{R}\). This looks pretty nasty! The equation is linear, but has non-constant coefficients and two parameters to account for. This equation is well understood analytically, but we’re only going to look at it computationally. In fact we can use some ideas from the third homework!

These functions arise in a few places throughout physics and engineering. Below we will cast this as an eigenvalue problem, similar to those you may find in quantum mechanics. Supposedly these functions also arise in control systems, though I have less experience there. If you would like to read more you can consult the wikipedia page, and the entry in the DLMF.

Re-casting the problem¶

When posed over the whole real line the above equation will have non-trivial solutions for all choices of \(a\) and \(q\). However, we will look for only the solutions that are \(2\pi-\) periodic. In this case we will choose \(q\) freely, and \(a\) will be fixed to a discrete (but still infinte) set of values. To make this more clear, we can re-write the above equation as an eigenvalue problem:

(2)¶\[\begin{split}\begin{cases} \left(2q\cos(2x) - \frac{d^2}{dx^2}\right)y = ay \\ y(0) = y(2\pi) \end{cases}\end{split}\]

where we can see that \(a\) and \(y\) form an eigenvalue/eigenfunction pair for the operator in the parentheses on the left.

How can we deal with that complicated looking differential operator? We already know that computers (mostly) can not do calculus, nor can they do infinitely many things. Let’s apply the same technique as in the second homework.

First, we’ll fix \(\mathbf{x}\) to be an array of equispaced grid points on \([0,2\pi)\). Then let \(\mathbf{T}\) be the discrete Fourier transformation matrix for this grid (as defined in HW3). A discrete version of the above equation can be written as:

(3)¶\[\left[\mathbf{V} - \mathbf{T}^{-1}\mathbf{D}_k\mathbf{T}\right]\mathbf{y} = a\mathbf{y}\]

where \(\mathbf{V}\) is a diagonal matrix with entries \(\lbrace2q\cos(2x_1),2q\cos(2x_2),2q\cos(2x_3),\ldots\rbrace\) and \(\mathbf{D}_k\) is a diagonal matrix with the squares of the wavenumbers as entries \(\lbrace 0,1,1,4,4,9,9,\ldots\rbrace\).

In the code we’ll pack everything inside the square braces into one matrix, \(\mathbf{H}\), and use LAPACK to obtain the eigenvalue decomposition.

The code¶

The code to solve this problem is split into several files:

mathieu.f90: This defines the main program and is responsible for calling the other modules to set things up, calling LAPACK, and writing the solutions to disc.
diffop.f90: This defines the diffop module which is responsible for building the matrices (discrete differential operators) we need.
problemsetup.f90: This holds the problem specific parameters, like the number of grid points, the value for \(q\), and what to label the output files. This also has routines to read an input file so that we don’t need to re-compile to try different parameter values.
utility.f90: This defines useful constants like \(pi\) and the precision.

There are also some helper files to make our life easier when using the code:

Makefile: This, of course, lets the program make do all of the compilation work for us.
mathieu.init: This defines the settings that get read in by the module from problemsetup.f90. Changing values here allow us to examine the effect of various parameter choices without needing to re-compile each time.
parStudy.sh: This is a very simple shell script that will batch out many jobs to perform a small parameter study.
plotEVals.py and plotEFuncs.py: These are basic python plotting scripts so that we can visualize the results from the code.

All of these can be obtained in this tarball which can be uncompressed by calling tar -xvzf MathieuFunctions.tar.gz.

Note

This example relies on a couple of the routines from HW3. You’ll need to insert your solutions from there for the code above to work!!

The `mathieu.f90` file¶

! File: mathieu.f90
! Author: Ian May
! Purpose: Compute Mathieu functions of the first kind numerically

program mathieu
  use utility, only : fp, pi
  use problemsetup, only : Npts, qIdx, outFile, problemsetup_Init
  use diffop, only : diffop_Mathieu
  implicit none
  
  ! Size of transform
  ! Problem data
  integer :: info,lwork                               ! dgeev return code and size of workspace
  real (fp) :: dx                                     ! grid spacing
  real (fp), allocatable :: x(:)                      ! grid points
  real (fp), allocatable :: wr(:), wi(:)              ! real and imaginary parts of e-vals
  real (fp), allocatable :: work(:)                   ! Workspace for dgeev
  real (fp), allocatable :: H(:,:), vr(:,:), vl(:,:)  ! operator, right e-vecs, left e-vecs

  ! Loop variables
  integer :: i,j
  
  ! Read in the input file
  call problemsetup_Init('mathieu.init')
  
  ! Set local variables
  lwork = 8*Npts
  dx = 2*pi/Npts
  
  ! Allocate arrays to match problem size
  allocate(x(Npts))
  allocate(wr(Npts))
  allocate(wi(Npts))
  allocate(work(lwork))
  allocate(H(Npts,Npts))
  allocate(vl(Npts,Npts))
  allocate(vr(Npts,Npts))
  
  ! Set up physical domain
  do i=1,Npts
    x(i) = (i-1)*dx
  end do
  
  ! Get Mathieu operator
  call diffop_Mathieu(qIdx,x,H)
  
  ! Use Lapack to find the eigenvalues/eigenvectors of H
  call dgeev('N','V',Npts,H,Npts,wr,wi,vl,1,vr,Npts,work,lwork,info)
  
  ! Open the output file and write the columns
  open(20,file=outFile,status="replace")
  do i=1,Npts
    write(20,*) x(i), wr(i), wi(i), (vr(i,j), j=1,Npts)
  end do
  close(20)
  
  ! Deallocate all used space
  deallocate(x)
  deallocate(wr)
  deallocate(wi)
  deallocate(work)
  deallocate(H)
  deallocate(vl)
  deallocate(vr)
end program mathieu

Download this code

The `diffop.f90` file¶

module diffop

  use utility, only : fp, pi
  
  implicit none
  private

  public :: diffop_Mathieu
  
contains

  ! subroutine: diffop_SecondDeriv
  ! purpose: Create second derivative operator over a given (periodic) grid
  ! inputs: x -- Rank 1 array holding grid positions, assumed periodic and equispaced
  ! outputs: D -- Square rank 2 array for the discrete second derivative operator
  subroutine diffop_SecondDeriv(x,D)
    implicit none
    real (fp), intent(in)     :: x(:)
    real (fp), intent(in out) :: D(:,:)
    
    ! Local variables
    integer :: M, N, i
    real (fp) :: om, dx, k
    real (fp) :: T(size(D,1),size(D,2)), Tinv(size(D,1),size(D,2))
    
    ! Set sizes and base wavenumber
    M=size(D,1)
    N=size(D,2)
    dx = x(2)-x(1)
    om = 2*pi/(N*dx)
    
    ! Fill transformation matrices
    ! T comes from HW3
    ! Tinv is already contracted with the diagonal scaling
    T(1,:) = ! Something?
    Tinv(:,1) = 0.0_fp
    do i=2,M,2
      k = i*om/2
      T(i,:) = ! Another ?
      Tinv(:,i) = -k**2*cos(k*x)
      if (i+1 <= M) then
        T(i+1,:) = ! One more...
        Tinv(:,i+1) = -k**2*sin(k*x)
      end if
    end do

    ! Last column of T needs re-scaling for even M
    if (mod(M,2)==0) then
      T(M,:) = T(M,:)/2.0_fp
    end if

    ! Multiply together to get spectral differentiation matrix
    ! Recall Tinv already has diagonal scaling in it
    D = matmul(Tinv,T)
  end subroutine diffop_SecondDeriv

  ! subroutine: diffop_Mathieu
  ! purpose: Create the Mathieu differential operator for a given index over a given
  !          (periodic) grid
  ! inputs:  q -- Mathieu index
  !          x -- Rank 1 array holding grid positions, assumed periodic and equispaced
  ! outputs: H -- Square rank 2 array for the discrete Mathieu operator
  subroutine diffop_Mathieu(q, x, H)
    implicit none
    real (fp), intent(in)     :: q
    real (fp), intent(in)     :: x(:)
    real (fp), intent(in out) :: H(:,:)
    
    ! Local variables
    integer :: N, j
    real (fp) :: V(size(x))
    N = size(x)
    
    ! Set H to the negative second derivative operator
    call diffop_SecondDeriv(x,H)
    H = -1.0_fp * H
    
    ! Add the Mathieu potential along the diagonal
    V = 2*q*cos(2*x)
    do j=1,N
      H(j,j) = H(j,j) + V(j)
    end do
  end subroutine diffop_Mathieu

end module diffop

Download this code

What are these allocate and deallocate statements all about? Since the size of the grid and matrices are specified inside an input file, the memory requirements for these arrays can’t be known at compile time. The allocate and deallocate commands dynamically make room in memory for these arrays at compile time. We’ll revisit this idea in the chapter on C.

The `problemsetup.f90` file¶

! File: problemsetup.f90
! Author: Ian May
! Purpose: Define problem specific information
! Comments: The read_initFile* subroutines were taken from Prof. Lee's
!           Newton's method example code

module problemsetup

  use utility, only : fp, maxFileLen, maxStrLen
  
  implicit none
  private
  
  integer, public :: Npts                              ! Number of grid points to use in the discretization
  real (fp), public :: qIdx                            ! Index of Mathieu functions we are looking for
  character(len=maxStrLen), public :: runName, outFile ! Name of the run and the output file the code should write

  public :: problemsetup_Init

contains

  ! subroutine: problemsetup_Init
  ! purpose: Set module variables to values defined in an input file
  ! inputs: inFile -- String with name of input file to be read
  ! outputs: <none>
  subroutine problemsetup_Init(inFile)
    implicit none
    character(len=*), intent(in) :: inFile
    
    ! Fill in default values
    Npts = 31
    qIdx = 36.0_fp
    
    ! Read problem resolution
    Npts = read_initFileInt(inFile,'num_points')
    
    ! Read the Mathieu index
    qIdx = read_initFileReal(inFile,'q_index')
    
    ! Set the name of the run from the number of points and index
    write(runName,"('Mathieu_',I0.3,'_',I0.2)") Npts,floor(qIdx)
    print *, 'Running problem: ', runName
    
    ! Set the output file, note that // does string concatenation
    outFile = 'data/' // trim(runName) // '.dat'
  end subroutine problemsetup_Init

  ! function: read_initFileReal
  ! purpose: Pull one real value from an input file
  ! inputs: inFile -- String holding the name of the input file
  !         varName -- String that names the variable, this must be first entry on a line
  ! outputs: varValue -- Real value that will hold the result from the input file
  function read_initFileReal(inFile,varName) result(varValue)

    implicit none
    character(len=*),intent(IN) :: inFile,varName
    real (fp) :: varValue

    integer :: i,openStatus,inputStatus
    real :: simInitVars
    character(len=maxStrLen) :: simCharVars
    integer :: pos1,pos2

    open(unit = 10, file=inFile, status='old',IOSTAT=openStatus,FORM='FORMATTED',ACTION='READ')

    do i=1,maxFileLen
      read(10, FMT = 100, IOSTAT=inputStatus) simCharVars
      pos1 = index(simCharVars,varName)
      pos2 = pos1 + len_trim(varName)
      if (pos2 > len_trim(varName)) then
        read(simCharVars(pos2 + 1:),*) simInitVars
        varValue = simInitVars
      endif
    end do

    close(10)

100 FORMAT(A, 1X, F3.1)

  end function read_initFileReal

  ! function: read_initFileInt
  ! purpose: Pull one integer value from an input file
  ! inputs: inFile -- String holding the name of the input file
  !         varName -- String that names the variable, this must be first entry on a line
  ! outputs: varValue -- Integer value that will hold the result from the input file
  function read_initFileInt(inFile,varName) result(varValue)

    implicit none
    character(len=*),intent(IN) :: inFile,varName
    integer :: varValue

    integer :: i,openStatus,inputStatus
    integer :: simInitVars
    character(len=maxStrLen) :: simCharVars
    integer :: pos1,pos2

    open(unit = 11, file=inFile, status='old',IOSTAT=openStatus,FORM='FORMATTED',ACTION='READ')

    do i=1,maxFileLen
      read(11, FMT = 101, IOSTAT=inputStatus) simCharVars
      pos1 = index(simCharVars,varName)
      pos2 = pos1+len_trim(varName)
      if (pos2 > len_trim(varName)) then
        read(simCharVars(pos2+1:),*)simInitVars
        varValue = simInitVars
      endif
    end do

    close(11)

101 FORMAT(A, 1X, I5)

  end function read_initFileInt

end module problemsetup

Download this code

The `utility.f90` file¶

! File: utility.f90
! Author: Ian May
! Purpose: Define useful constants

module utility
  
  implicit none
  
  integer, parameter :: fp = selected_real_kind(15)
  integer, parameter :: maxFileLen=50, maxStrLen=200
  real (fp), parameter :: pi = acos(-1.0_fp)

contains
  
end module utility

Download this code

The `Makefile`¶

FC = gfortran
FFLAGS = -Wall -Wextra -Wno-surprising -fPIC -fmax-errors=1 -g -fcheck=all -fbacktrace -J ./obj
LFLAGS = -llapack

ODIR=obj

_OBJ = utility.o problemsetup.o diffop.o mathieu.o
OBJ = $(patsubst %,$(ODIR)/%,$(_OBJ))

mathieu.ex: $(OBJ)
	$(FC) -o $@ $^ $(LFLAGS) $(FFLAGS)

$(ODIR)/%.o: %.f90
	@mkdir -p $(ODIR)
	$(FC) -c -o $@ $< $(FFLAGS)

.PHONY: clean

clean:
	rm -f mathieu.ex obj/* *~

Download this code

Note the special flag -J ./obj in the Fortran flags section. This redirects all object and module files into the obj/ directory to keep our build directory somewhat clean. You can try putting the *.f90 files into a src/ subdirectory as well.

Fortran Example – Mathieu functions and characteristics¶

Re-casting the problem¶

The code¶

The `mathieu.f90` file¶

The `diffop.f90` file¶

The `problemsetup.f90` file¶

The `utility.f90` file¶

The `Makefile`¶

Table of Contents

Previous topic

Next topic

This Page

Fortran Example – Mathieu functions and characteristics¶

Re-casting the problem¶

The code¶

The mathieu.f90 file¶

The diffop.f90 file¶

The problemsetup.f90 file¶

The utility.f90 file¶

The Makefile¶

The `mathieu.f90` file¶

The `diffop.f90` file¶

The `problemsetup.f90` file¶

The `utility.f90` file¶

The `Makefile`¶