The discussion so far has been for first order ordinary differential equations. All the methods given may be applied to higher ordinary differential equations, provided it is possible to write an explicit expression for the highest order derivative and the system has a complete set of initial conditions. Consider some equation

(96)

,

where at t = t0 we know the values of y, dy/dt, d2y/dt2, …, dn1y/dtn1. By writing

x0=y, x1=dy/dt, x2=d2y/dt2, …, xn-1=dn-1y/dtn-1, we may express this as the system of equations

x0‘ = x1

x1‘ = x2x2‘ = x3

. . . .

x ‘ = x

n-2 n-1

xn-1‘ = f(t,x0,x1,…,xn-2), (97) and use the

standard methods for updating each xi for some tn+1 before

proceeding to the next time st p. A decision needs to be made as to whether the values of xi for tn or tn+1 are o be used on the right hand side of the

equation for xn-1‘. This decision may affec the order and convergence of the method. Detailed analysis may be undertaken in a manner similar to that for the first order ordinary differential quations.

7.2 Boundary value problems

For second (and higher) order odes, two (or more) initial/boundary conditions are required. If these two conditions do not correspond to the same point in time/space, then the simple extension of the first order methods outlined in section 7.1 can not be applied without modification. There are two relatively simple approaches to solve such equations.

7.2.1 Shooting method

Suppose we are solving a second order equation of the form y” = f(t,y,y’) subject to y(0) = c0 and y(1) = c1. With the shooting method we apply the y(0)=c0 boundary condition and make some guess that y’(0) = α0. This gives us two initial conditions so that we may apply the simple time-steppingmethods already discussed in section7.1 . The calculation proceeds until we have a value

y” + ay‘ + by = c,		(98)
with boundary conditions y(t0) = α and y’(t1) = β. If we use the central
difference approximations
y’i =~ (Yi+1 – Yi-1)/2 t,		(99a)
y”i =~ (Yi+1 – 2Yi + Yi-1)/ t2,		(	99	b)
we can write the system as
Y0 = α,
(1+1/2a t)Y0 + (b t2–2)Y1 + (1-1/2a	t)Y2 = c	t2,
(1+1/2a t)Y1 + (b t2–2)Y2 + (1-1/2a t)Y3 = c t2,
…
(1+1/2a t)Yn-1 + (b t2–2)Yn-1 + (1-1/2a	t)Yn = c	t2,
Yn – Yn-1 = βΔt		(100)

Newton-Raphson algorithm, similar in principle to that introduced in section 3.4 but where a system of linear equations requires solution for each iteration.

system will have [(m+1)(n+1)]2 elements, the storage and computational cost of such a solution will become prohibitive even for relatively modest m and n.

The structure of the system ensures A is relatively sparse, consisting of a tridiagonal core with one nonzero diagonal above and another below this. These nonzero diagonals are offset by either m or n from the leading diagonal. Provided pivoting (if required) is conducted in such a way that it does not place any nonzero elements outside this band then solution by Gauss Elimination or LU Decomposition will only produce nonzero elements inside this band, substantially reducing the storage and computational requirements (see section4.4 ). Careful choice of the order of the matrix elements (i.e. by x or by

y) may help reduce the size of this matrix so that it need contain only O(m3) elements for a square domain.

Because of the wide spread need to solve Laplace’s and related equations, specialised solvers have been developed for this problem. One of the best of these is Hockney’s method for solving Ax = b which may be used to reduce a block tridiagonal matrix (and the corresponding right-hand side) of the form

(10

4)

,

into a block diagonal matrix of the form

(10

5)

,

where and are themselves block tridiagonal matrices and I is an identiy matrix.. This process may be performed iteratively to reduce an n dimensional finite difference approximation to Laplace’s equation to a tridiagonal system of equations with n-1 applications. The computational cost is O(p log p), where p is

the total number of mesh points. The main drawback of this method is that the boundary conditions must be able to be cast into the block tridiagonal format.

8.1.2 Relaxation

An alternative to direct solution of the finite difference equations is an iterative numerical solution. These iterative methods are often referred to as relaxation methods as an initial guess at the solution is allowed to slowly relax towards the true solution, reducing the errors as it does so. There are a variety of approaches with differing complexity and speed. We shall introduce these methods before looking at the basic mathematics behind them.

8.1.2.1 Jacobi

The Jacobi Iteration is the simplest approach. For clarity we consider the special case when x = y. To find the solution for a two-dimensional Laplace equation simply:

1.Initialise Φij to some initial guess.

1.Replace old solution Φ with new estimate Φ*.

The coefficients in the expression (here all 1/4) used to calculate the refined estimate is often ref rr d t as the stencil or template. Higher order approximations may b obtained by simply employing a stencil which utilises more points. Other quations (e.g. the bi-harmonic equation, 4Ψ = 0) may be solved by introduc ng a stencil appropriate to that equation.

While very simple and cheap per iteration, the Jacobi Iteration is very slow to converge, especially for larger grids. Corrections to errors in the estimate Φijdiffuse only slowly from the boundaries taking O(max(m,n)) iterations to diffuse across the entire mesh.

8.1.2.2 Gauss-Seidel

The Gauss-Seidel Iteration is very similar to the Jacobi Iteration, the only difference being that the new estimate Φ*ij is returned to the solution Φij as soon as it is completed, allowing it to be used immediately rather than deferring its use to the next iteration. The advantages of this are:

●Less memory required (there is no need to store Φ*).

On the other hand, the method is less amenable to vectorisation as, for a given iteration, the new estimate of one mesh point is dependent on the new estimates for those already scanned.

8.1.2.3 Red-Black ordering

A variant on the Gauss-Seidel Iteration is obtained by updating the solution Φijin two passes rather than one. If we consider the mesh points as a chess board, then the white squares would be updated on the first pass and the black squares on the second pass. The advantages

It has been found that the errors in the solution obtained by any of the three preceding methods decrease only slowly an often decrease in a monotonic manner. Hence, rather than setting

for some value σ. Th optimal value of σ will depend on the problem being solved and may vary as the iteration process converges. Typically, however, a value of around 1.2 to 1.4 produces good results. In some special cases it is possible to determ ne an optimal value analytically.

8.1.3 Multigrid*

The big problem with relaxation methods is their slow convergence. If σ = 1 then

application of the stencil removes all the error in the solution at the wave length

of the mesh for that point, but has little impact on larger wave lengths. This may

be seen if we consider the one-dimensional equation d2/dx2 = 0 subject to

ϕ(x=0) = 0 and ϕ(x=1) = 1. Suppose our initial guess for the iterative solution is

that Φi = 0 for all internal mesh points. With the Jacobi Iteration the correction to

the internal points diffuses only slowly along from x = 1.

Multigrid methods try to improve the rate of convergence by considering the problem of a hierarchy of grids. The larg r wave length errors in the solution are dissipated on a coarser grid while the shorter wave length errors are dissipated on a finer grid. for the example con id r d above, the solution would converge in one complete Jacobi multigrid iteration, compared with the slow asymptotic convergence above.

For linear problems, the basic multigrid algorithm for one complete iteration may be described as

1.Select the initial finest grid resolution p=P0 and set b(p) = 0 and make some initial guess at the solution Φ(p)

2.If at coarsest resolution (p=0) then solve A(p)Φ(p)=b(p) exactly and jump to step 7

5.Coarsen the error b(p-1)r to the next coarser grid and decrement p

7.Refine the correction to the next finest grid Φ(p+1) = Φ(p+1)+αΦ(p)and increment p

9.If not at current finest grid (P0), repeat from step 7

10.If not at final desired grid, increment P0 and repeat from step 7

Typically the relaxation steps will be performed using Successive Over Relaxtion with Red-Black ordering and some relaxation coefficient σ. The hierarchy of grids is normally chosen to differ in dimensions by a factor of 2 in each direction. The factor α is typically less than unity and effectively damps possible instabilities in the convergence. The refining of the correction to a finer grid will be achieved by (bi-)linear or higher order interpolation, and the coarsening may simply be by sub-sampling or averaging the error vector r.

It has been found that the number of iterations required to reach a given level of convergence is more or less independent of the number of mesh points. As the number of operations per complete iteration for n mesh points is O(n)+O(n/2d)+ +O(n/22d)+…, where d is the number of dimensions in the problem, then it can be seen that the Multigrid method may often be faster than a direction solution (which will require O(n3), O(n2) or O(n log n) operations, depending on the method used). This is particularly true if n is large or there are a large number of dimensions in the problem. For small problems, the coefficient in front of the n for the Multigrid solution may be relatively large so that direct solution may be faster.

A further advantage of Multigrid and oth r iterative methods when compared with direct solution, is that irregular shaped domains or complex boundary conditions are implemented more easily. The difficulty with this for the Multigrid method is that care must be taken in ord r to ensure consistent boundary conditions in the embedded problems.

8.1.4 The mathematics of relaxation*

In principle, relaxation methods which are the basis of the Jacobi, Gauss-Seidel,Successive Over R laxation and Multigrid methods may be applied to any system of linear equations to interatively improve an approximation to the exact solution. The bas s for this is identical to the Direct Iteration method described in section3.6 . We start by writing the vector function

f(x) = Ax b,

and search for the vector of roots to f(x) = 0 by writing

xn+1 = g(xn), 9)

with D a diagonal matrix (zero for all off-diagonal elements) which may be chosen arbitrarily. We may analyse this system by following our earlier analysis for the Direct Iteration method (section3.6 ). Let us assume the exact solution is x* = g(x*), then

ε = x x*

n+1 n+1

=D1{[A+D]xn b} D1{[A+D]x* b}

=D1[A+D](xn x*)

=D1[A+D]εn

= {D1[A+D]}n+1 ε0.

From this it is clear that convergence will be linear and requires

where B = D1[A+D] for some suitable norm. As any error vector εn may be written as a linear combination of the eigen vectors of our matrix B, it is sufficient for us to consider the eigen value problem

(11

Bεn = λεn,2)

and require max(|λ|) to be less than unity. In the asymptotic limit, the smaller the magnitude of this maximum eigen valu the more rapid the convergence. The convergence remains, however, linear

Since we have the ability to choo e the diagonal matrix D, and since it is the eigen values of B = D1[A+D] ra h r than A itself which are important, careful choice of D can aid the speed a which the method converges. Typically this means selecting D so that the diagonal of B is small.

8.1.4.1 Jacobi a d Gauss-Seidel for Laplace equation*

The structure of the finite difference approximation to Laplace’s equation lends itself to these relaxation methods. In one dimension,

(11

3)

and both Jacobi and Gauss-Seidel iterations take D as 2I (I is the identity matrix) on the diagonal to give B = D1[A+D] as

(11

4)

The eigen values λ of this matrix are given by the roots of

It can be shown that for a system of any size following this general form, all the eigen values satisfy |λ| < 1, thus proving the relaxation method will always converge. As we increase the number of mesh points, the number of eigen values increases and gradually fills up the range |λ| < 1, with the numerically largest eigen values becoming closer to unity. As a result of λ1, the convergence of the relaxation method slows considerably for large problems. A similar analysis may be applied to Laplace’s equation in two or more dimensions, although the expressions for the determinant and eigen values is correspondingly more complex.

The large eigen values are responsible for decreasing the error over large distances (many mesh points). The multigrid approach enables the solution to

converge using a much smaller system of equations and hence smaller eigen values for the larger distances, bypassing the slow convergence of the basic relaxation method.

8.1.4.2 Successive Over Relaxation for Laplace equation*

The analysis of the Jacobi and Gauss-Seidel iterations may be applied equally well to Successive Over Relaxation. The main difference is that D = (2/σ)I so that

(119)

and the corresponding eigen values are related by (1σλ)2 equal to the values tabulated above. Thus if σ is chosen inappropriately, the eigen values of B will exceed unity and the relaxation method will iverge. On the otherhand, careful choise of σ will allow the eigen values of B to be less than those for Jacobi and Gauss-Seidel, thus increasing the rat of convergence.

8.1.4.3 Other equations*

Relaxation methods may be appli d to other differential equations or more general systems of linear equa ions in a similar manner. As a rule of thumb, the solution will converge if the A matrix is diagonally dominant, i.e. the numerically largest values occur o the diagonal. If this is not the case, SOR can still be used, but it may be n cessary to choose σ < 1 whereas for Laplace’s equation

σ>= 1 produces a better rate of convergence.

8.1.5 FFT*

One of the most common ways of solving Laplace’s equation is to take the Fourier transform of the equation to convert it into wave number space and there solve the resulting algebraic equations. This conversion process can be very efficient if the Fast Fourier Transform algorithm is used, allowing a solution to be evaluated with O(n log n) operations.

In its simplest form the FFT algorithm requires there to be n = 2p mesh points in the direction(s) to be transformed. The efficiency of the algorithm is achieved by first calculating the transform of pairs of points, then of pairs of transforms, then

of pairs of pairs and so on up to the full resolution. The idea is to divide and conquer! Details of the FFT algorithm may be found in any standard text.

The Poisson equation 2φ = f(x) may be treated using the same techniques as Laplace’s equation. It is simply necessary to set the right-hand side to f, scaled suitably to reflect any scaling in A.

8.3 Diffusion equation

Consider the two-dimensional diffusion equation,

(12

subject to u(x,y,t) = 0 on the boundaries x=0,1 and y=0,1. Suppose the initial conditions are u(x,y,t=0) = u0(x,y) and we wish to evaluate the solution for t > 0. We shall explore some of the options for achieving this in the following sections.

8.3.1 Semi-discretisation

One of the simplest and mos us ful approaches is to discretise the equation in space and then solve a system of (coupled) ordinary differential equations in time in order to calculate the solution. Using a square mesh of step size

x = y = 1/m, and taki g the diffusivity D = 1, we may utilise our earlier approximation f r th Laplacian operator (equation ( 103 )) to obtain

(12

1) for the internal points i=1,m-1 and j=1,m-1. On the boundaries (i=0,j), (i=m,j), (i,j=0) and (i,j=m) we simply have uij=0. If Uij represents our approximation of u at the mesh points xij, then we must simply solve the (m-1)2 coupled ordinary differential equations

In principle we may utilise any of the time stepping algorithms discussed in earlier lectures to solve this system. As we shall see, however, care needs to be taken to ensure the method chosen produces a stable solution.

			(12
	where the Courant number	,	3)
	μ =	t/ x2,	(12
			4)

U(n)i,j = sin(αi)sin(βj) {1 4μ[sin2(α/2) + sin2(β/2)]}n,

(127)

which then requires |1 4μ[sin2(α/2) + sin2(β/2)]| < 1 for this to converge as n>infinity. For this to be satisfied for arbitrary α and β we require μ < 1/4. Thus we must ensure

A doubling of the spatial resolution therefore requires a factor of four more time steps so overall the expense of the computation increases sixteen-fold.

The analysis for the diffusion equation in one or three dimensions may be computed in a similar manner.

8.3.4 Model for general initial conditions