INSTRUCTIONS FOR SOFTWARE USE Part IV, problem (2.2) ============================================================== The file matrixroutines.f contains four fortran-77 linear algebra subroutines that you may find useful to set up an inverse iteration code according to the pseudo algorithm listed in the class notes (part IV, p. XXX). The routines are: tribu --> triangularizes a column-shifted matrix (see section 1 below) rhsbu --> uses triangularized form of matrix to solve for a vector of unknown, given a known RHS vector fixbc --> implements Dirichlet-nul boundary conditions at endpoints (see section 2 below) matxvec --> performs a matrix-vector multiply Those of you who will bother to look into the routines will notice that they are more complicated than they could have been; this is because they are in fact routines for general banded systems, being specialized to the tridiagonal case with zero-endpoint boundary conditions. You may want to hang on to these routines even ater the class is over; they may end up useful for other problems you may encounter in the future. ============================================================== 1. Storing the matrices Discretization of a PDE using finite differences in 1-D typically produces BANDED matrices. For example, using the second-order centered finite differences given in equations (2.XX) on a mesh of N gridpoints will produce a tridiagonal NxX matrix: [ D U ] [ L D U ] [ L D U ] [ L D U ] [ . . . ] [ L D U ] [ L D ] where all unmarked entries are zeros. For example, a simple second derivative on a constant-increment 1-D mesh would yield L=U=1/(Dx)**2 and D=-2/(Dx)**2, except for the first and last node of the mesh (mapping onto the first and last line in the matrix, but we'll deal with that later). Clearly it doesn't make much sense to store this as a NxN array when we only have 3xN-2 non-zero elements. This matrix has a FULL BANDWIDTH (mb) of 3, and a HALF BANDWIDTH (mhb) of 2. Consider the shifting operation produced my moving the matrix elements (i,j) vertically according to the shifting rule (i,j) --> (i+mhb-j,j) This produces the shifted version of [ D U ] [ U U U U . U ] [ L D U ] [ D D D D . D D ] [ L D U ] [ L L L . L L L ] [ L D U ] --> [ ] [ . . . ] [ ] [ L D U ] [ ] [ L D ] [ ] which is conveniently stored as a (mhb x N) array The price to pay for this saving in storage is that matrix operations, even simple ones such as matrix-vector multiply, done using the shifted version of the matrix, will need to do some additional integer algebra to keep track of the real position of each element in the original tridiagonal square matrix. ============================================================== 2. Hardwiring the boundary conditions Consider the matrix system [ D U ]{a1} {f1} [ L D U ]{a2} {f2} [ L D U ]{a3} {f3} [ L D U ]{a4} = {f4} [ . . . ]{. } {. } [ L D U ]{. } {. } [ L D ]{aN} {fN} where the {a} vector is to be solved for, given the known RHS vector {f}. Suppose that we need to impose the boundary conditions a1=A aN=B on this system of equationss. This can be done by the following modifications to the matrix system: [ 1 0 ]{a1} {A } [ 0 D U ]{a2} {f2 -L A} [ L D U ]{a3} {f3 } [ L D U ]{a4} = {f4 } [ . . . ]{. } {. } [ L D 0 ]{. } {fN-1 -U B} [ 0 1 ]{aN} {B } Clearly, inversion of this sysmtem will yield a1=A and aN=B. If the boundary consitions are a1=aN=0, which is the case for this problem, then this becomes simply [ 1 0 ]{a1} {0 } [ 0 D U ]{a2} {f2} [ L D U ]{a3} {f3} [ L D U ]{a4} = {f4} [ . . . ]{. } {. } [ L D 0 ]{. } {. } [ 0 1 ]{aN} {0 } Once again because the matrix is stored in the form of rectangular array, the zeroing of the first and last lines and columns involves some integer algebra to take into account the shifting. Subroutine FIXBC carries out these modifications for you, for the special case of A=B=0. ============================================================== 3. Solving a system of algebraic equations Consider now the algebraic system defined by [G]{a}={f} where [G] is the (mhb x N) representation of a tridiagonal (N x N) square matrix {f} is a vector (dimension N) of known quantities {a} is a vector (dimension N) of unknown to be solved for subject to the boundary conditions a(1)=0 and a(N)=0; this is the sort of thing you have on your hands in using inverse iteration to to the toroidal diffusive eigenmode problem. With the matriz and vector sizes properly defined, the following sequences of calls will solve this system: mhb=2 mb=3 call fixbc(N,mhb,mb,G) call tribu(N,mhb,mb,G) call rhsbu(N,mhb,mb,G,f,a) And that's it! Note that the subroutines are not restricted to tridiagonal systems (mhb=3, mb=2); this is why you must input explicitly the full- and half-bandwith values. ============================================================== 4. Some hints for coding up inverse iteration You may recall from the discussion of section 2.1.2 that the inverse iteration algorithm requires a reasonable guess at the sought-after eigenvalue, as well as (less critically) the sought-after eigenvector. For guesses at the eigenvalues you may want to use the eigenvalues obtained for poloidal fields (Figure 2.1) For guesses at eigenvectors, you may want to try something like sin( (n+1) * pi * r) 0