.. _ch02-fortran-mathieu: ========================================================= 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, :math:`y(x)`, are defined as the solutions of the differential equation: .. math:: \frac{d^2y}{dx^2} + \left(a-2q\cos(2x)\right)y(x) = 0 where :math:`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 :math:`a` and :math:`q`. However, we will look for only the solutions that are :math:`2\pi-` periodic. In this case we will choose :math:`q` freely, and :math:`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: .. math:: \begin{cases} \left(2q\cos(2x) - \frac{d^2}{dx^2}\right)y = ay \\ y(0) = y(2\pi) \end{cases} where we can see that :math:`a` and :math:`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 :math:`\mathbf{x}` to be an array of equispaced grid points on :math:`[0,2\pi)`. Then let :math:`\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: .. math:: \left[\mathbf{V} - \mathbf{T}^{-1}\mathbf{D}_k\mathbf{T}\right]\mathbf{y} = a\mathbf{y} where :math:`\mathbf{V}` is a diagonal matrix with entries :math:`\lbrace2q\cos(2x_1),2q\cos(2x_2),2q\cos(2x_3),\ldots\rbrace` and :math:`\mathbf{D}_k` is a diagonal matrix with the squares of the wavenumbers as entries :math:`\lbrace 0,1,1,4,4,9,9,\ldots\rbrace`. In the code we'll pack everything inside the square braces into one matrix, :math:`\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 :math:`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 :math:`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 :download:`this tarball <./codes/MathieuFunctions.tar.gz>` 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 ------------------------- .. literalinclude:: ./codes/MathieuFunctions/mathieu.f90 :language: fortran :linenos: :download:`Download this code <./codes/MathieuFunctions/mathieu.f90>` The ``diffop.f90`` file ------------------------ .. literalinclude:: ./codes/MathieuFunctions/diffop.f90 :language: fortran :linenos: :download:`Download this code <./codes/MathieuFunctions/diffop.f90>` 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 ------------------------------ .. literalinclude:: ./codes/MathieuFunctions/problemsetup.f90 :language: fortran :linenos: :download:`Download this code <./codes/MathieuFunctions/problemsetup.f90>` The ``utility.f90`` file ------------------------- .. literalinclude:: ./codes/MathieuFunctions/utility.f90 :language: fortran :linenos: :download:`Download this code <./codes/MathieuFunctions/utility.f90>` The ``Makefile`` ------------------ .. literalinclude:: ./codes/MathieuFunctions/Makefile :language: make :linenos: :download:`Download this code <./codes/MathieuFunctions/Makefile>` 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.