7. numerical method Higher order ordinary differential equations numerical method
 Posted by Engineering Helpline Pvt. Ltd
 Categories Fifth Semester
 Date March 31, 2018
 Comments 0 comment
7.Higher order ordinary differential equations
7.1Initial value problems
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, …, xn1=dn1y/dtn1, we may express this as the system of equations
x0‘ = x1
x1‘ = x2x2‘ = x3
. . . .
x ‘ = x
n2 n1
xn1‘ = f(t,x0,x1,…,xn2), (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 xn1‘. 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 timesteppingmethods already discussed in section7.1 . The calculation proceeds until we have a value
for y(1). If this does not satisfy y(1) = c1 to some acceptable tolerance, we revise our guess for y’(0) to some value α1, say, and repeat the time integration to obtain an new value for y(1). This process continues until we hity(1)=c1 to the acceptable tolerance. The number of iterations which will need to be made in order to achieve an acceptable tolerance will depend on how good the refinement algorithm for α is. We may use the root finding methods discussed in section 3 to undertake this refinement.
The same approach can be applied to higher order ordinary differential equations. For a system of order n with m boundary conditions at t = t0 and
nm boundary conditions at t = t1, we will require guesses for nm initial conditions. The computational cost of refining these nm guesses will rapidly become large as the dimensions of the space increase.
7.2.2 Linear equations
The alternative is to rewrite the equations using a finite difference approximation with step size t = (t1–t0)/N to produce a system of N+1 simultaneous equations. Consider the second order linear system
y” + ay‘ + by = c,

(98)


with boundary conditions y(t0) = α and y’(t1) = β. If we use the central


difference approximations


y’i =~ (Yi+1 – Yi1)/2 t,

(99a)


y”i =~ (Yi+1 – 2Yi + Yi1)/ t2,

(

b)


we can write the system as


Y0 = α,


(1+1/2a t)Y0 + (b t2–2)Y1 + (11/2a

t)Y2 = c

t2,


(1+1/2a t)Y1 + (b t2–2)Y2 + (11/2a t)Y3 = c t2,


…


(1+1/2a t)Yn1 + (b t2–2)Yn1 + (11/2a

t)Yn = c

t2,


Yn – Yn1 = βΔt

(100)

This tridiagonal system may be readily solved using the method
discussed in section4.5 .
Higher order linear equations may be catered for in a similar manner and the matrix representing the system of equations will remain banded, but not as sparse as tridiagonal. The solution may be undertaken using the modified LU decomposition introduced in section4.4 .
Nonlinear equations may also be solved using this approach, but will require an iterative solution of the resulting matrix system Ax = b as the matrix A will be a function of x. In most circumstances this is most efficiently achieved through a
NewtonRaphson algorithm, similar in principle to that introduced in section 3.4 but where a system of linear equations requires solution for each iteration.
7.3Other considerations*
7.3.1Truncation error*
7.3.2Error and step control*
8.Partial differential equations
8.1Laplace equation
Consider the Laplace equation in two dimensions
(101)
,
in some rectangular domain described by x in [x0,x1], y in [y0,y1]. Suppose we discretise the solution onto a m+1 by n+1 rectangular grid (or mesh) given by
xi = x0 + i x, yj = y0 + j y where i=0,m, j=0,n. The mesh spacing is
x= (x1x0)/m and y = (y1y0)/n. Let ij = (xi,yj) be the exact solution at the mesh point i,j, and Φij =~ ij be the approximate solution at that mesh point.
By considering the Taylor Series expansion for about some mesh point i,j,
(102a)
,
,
( 102
b)
it is clear that we may appr ximate 2ϕ/x2 and 2ϕ/y2to the first order using the four adjacent mesh points to obtain the finite difference approximation
(103)
for the internal points 0<i<m, 0<j<n. In addition to this we will have either Dirichlet, von Neumann or mixed boundary conditions to specify the boundary values of ij. The system of linear equations described by ( 103 ) in combination
with the boundary conditions may be solved in a variety of ways.
8.1.1 Direct solution
Provided the boundary conditions are linear in , our finite difference approximation is itself linear and the resulting system of equations may be solved directly using Gauss Elimination as discussed in section4.1 . This approach may be feasible if the total number of mesh points (m+1)(n+1) required is relatively small, but as the matrix A used to represent the complete
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 righthand 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 n1 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 twodimensional Laplace equation simply:
1.Initialise Φij to some initial guess.
2.Apply the boundary conditions.
3.For each internal mesh point set
Φ* = (Φ

+ Φ

+ Φ

+ Φ )/4.

(10


6)


ij

i+1,j

i1,j i,j+1

i,j1

1.Replace old solution Φ with new estimate Φ*.
2.If solution does not sa isfy olerance, repeat from step 2.
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 biharmonic 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 GaussSeidel
The GaussSeidel 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 Φ*).
●Faster convergence (although still relatively slow).
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 RedBlack ordering
A variant on the GaussSeidel 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
●No interdependence of the solution updates within a single pass aids vectorisation.
●Faster convergence at low wave numbers.
8.1.2.4Successive Over Relaxation (SOR)
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
Φ* = (Φ

+ Φ

+ Φ

+ Φ )/4,


ij

i+1,j

i1,j

i,j+1

i,j1


for each internal mesh point, w u e


Φ* = (1σ)Φ + σ(Φ + Φ

+ Φ

(10


ij

ij

i+1,j

i1,j

i,j+1 + Φi,j1)/4,

7)


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 onedimensional 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
3.Relax the solution at the current grid resolution, applying boundary conditions
4.Calculate the error r = AΦ(p)b(p)
5.Coarsen the error b(p1)r to the next coarser grid and decrement p
6.Repeat from step 2
7.Refine the correction to the next finest grid Φ(p+1) = Φ(p+1)+αΦ(p)and increment p
8.Relax the solution at the current grid resolution, applying boundary conditions
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
11.If not converged, repeat from step 2.
Typically the relaxation steps will be performed using Successive Over Relaxtion with RedBlack 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 subsampling 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, GaussSeidel,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)
where


g(x) = D1{[A+D]x b},

(11


0)


with D a diagonal matrix (zero for all offdiagonal 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
εn+1 = Bεn < εn,

(111)

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 GaussSeidel 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 GaussSeidel 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
det(BλI) = 0.

(11


5)


In this case the determinant may be obtained using the recurrence relation


det(Bλ)(n) = λ det(Bλ)(n1)

1/4 det(Bλ)(n2)

,

(11


6)


where the subscript gives the size of the matrix B. From this we may see


det(Bλ)(1) = λ ,


det(Bλ)(2) = λ2

¼ ,


det(Bλ)(3) = λ3 + ½λ ,


det(Bλ)(4) = λ4 ¾ λ2 + (1/16) ,


det(Bλ)(5) = λ5 + λ3

(3/16)λ ,


det(Bλ)(6) = λ6 (5/4)λ4 + (3/8)λ2 (1/64) ,


…

(117)


which may be s lv d to give the eigen values


λ(1) = 0 ,


λ2(2) = 1/4 ,


λ2(3) = 0, 1/2 ,


λ2(4) = (3 5)/8 ,


λ2(5) = 0, 1/4, 3/4 ,


…

(118)

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 GaussSeidel 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 GaussSeidel, 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.
8.1.6Boundary elements*
8.1.7Finite elements*
8.2Poisson equation
The Poisson equation 2φ = f(x) may be treated using the same techniques as Laplace’s equation. It is simply necessary to set the righthand side to f, scaled suitably to reflect any scaling in A.
8.3 Diffusion equation
Consider the twodimensional diffusion equation,
(12
,

0)

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 Semidiscretisation
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,m1 and j=1,m1. 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 (m1)2 coupled ordinary differential equations
(12


.

2)

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.
8.3.2 Euler method
Applying the Euler method Yn+1 = Yn+ tf(Yn,tn) to our spatially discretised diffusion equation gives
(12


where the Courant number

,

3)


μ =

t/ x2,

(12


4)


describes the size of the time step relative to the spatial discretisation. As we shall see, stability of the solution depends on μ in contrast to an ordinary differential equation where it is a function of the time step t only.
8.3.3 Stability
Stability of the Euler method solving the diffusio equation may be analysed in a similar way to that for ordinary differential eq atio s. We start by asking the question “does the Euler method converge as t>i finity?” The exact solution will have u > 0 and the numerical solution m st also do this if it is to be stable.
We choose
U(0) i,j=sin(αi) sin(βj),
for some α and β chosen as mul iples of π/m to satisfy u = 0 on the boundaries. Substituting this into ( 123 ) gives
U(1) i,j = sin(αi)sin(βj) + μ{sin[α(i+1)]sin(βj) + sin[α(i1)]sin(βj)
+sin(α )sin[β(j+1)] + sin(αi)sin[β(j1)] 4 sin(αi)sin(βj)}
=sin(αi)sin(βj) + μ{[sin(αi)cos(α) + cos(αi)sin(α)]sin(βj) + [sin(αi)cos(α) cos(αi)sin(α)]sin(βj)
+sin(αi)[sin(βj)cos(β) + cos(βj)sin(β)] + sin(αi)[sin(βj)cos(β) cos(βj)sin(β)] 4 sin(αi)sin(βj)}
=sin(αi)sin(βj) + 2μ{sin(αi)cos(α) sin(βj) + sin(αi)sin(βj)cos(β) 2 sin(αi)sin(βj)}
=sin(αi)sin(βj){1 + 2μ[cos(α) + cos(β) 2]}
=sin(αi)sin(βj){1 4μ[sin2(α/2) + sin2(β/2)]}.
Applying this at consecutive times shows the solution at time tn is
(12
5)
(126)
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
t < x2/4.

(12


8)


A doubling of the spatial resolution therefore requires a factor of four more time steps so overall the expense of the computation increases sixteenfold.
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
You may also like
eartquake protection:building technology
1 June, 2019
numerical method note
9 February, 2019
chapter wise note concrete technology
22 January, 2019