Finite Difference Method (FDM) Chaiwoot Boonyasiriwat August 27, - - PowerPoint PPT Presentation

Finite Difference Method (FDM)

Chaiwoot Boonyasiriwat

August 27, 2020

2

▪ The value of a derivative of function can be approximated by finite difference formulas based on function values at discrete points. ▪ Differential equations can then be solved by replacing derivatives with finite difference approximations. ▪ This discretization is called finite difference method. ▪ This leads to a system of algebraic equations which can be solved using numerical methods on a computer. ▪ A numerical solution from FDM are only known at discrete points in space and/or time. ▪ Accuracy of the solution depends on the order of accuracy of the finite difference approximations used.

Introduction to FDM

▪ Let u(x) be a function of one variable. ▪ A derivative of u(x) can be approximated based on the values of the function at discrete points. ▪ Recall the definition of derivative ▪ If the limit is dropped out, we then obtain a formula for approximating the first derivative ▪ This formula is a one-sided finite difference approximation.

Finite Difference Approximation

▪ Forward finite-difference scheme ▪ Backward finite-difference scheme ▪ Centered finite-difference scheme

Commonly Used Formulas

Example 1

LeVeque (2007, p. 4)

D0u(x) is obviously more accurate than D-u(x) and D+u(x)

Another finite-difference scheme

Asymmetric Approximation

First-order accurate approximation First-order accurate approximation Second-order accurate approximation Third-order accurate approximation

Order of Accuracy

Let u(x) = sin(x). Approximating u'(x) at x = 1 using various values of step length h, we obtain numerical errors from each approximation compared to the exact solution:

Example 2

LeVeque (2007, p. 4-5)

According to the results, we can see that In general, the approximation error can be written as where p is the order of approximation.

Example 2 (continued)

LeVeque (2007, p. 5)

Absolute errors versus h on a log-log scale

Example 2 (continued)

LeVeque (2007, p. 6)

h |E|

Taylor’s series expansion is the standard approach for analyzing the truncation error of a finite-difference approximation.

Truncation Error

12

Rearranging the first expansion yields

Truncation Error (continued)

Truncation error

13

Big O Notation: Let f and g be real-valued functions. One can write if and only if there exists a positive real number M and a real number x0 such that If a truncation error is then

Big O Notation

https://en.wikipedia.org/wiki/Big_O_notation

14

To make finite difference expressions concise and simple to handle, we usually use the following index notations.

Index Notation and Convention

15

▪ FD stencil is a graphical representation of a finite difference approximation. ▪ Stencil makes a complicated FD expression looks simpler.

FD Stencil

i - 2 i - 1 i i + 1 i + 2

• 1/12 4/3 -5/2 4/3 -1/12

16

▪ Finite difference coefficients can be derived using the method of undetermined coefficients. ▪ Example 3: Suppose we want to derive a one-sided FD approximation to u'(x) based on the function values u(x), u(x – h), and u(x – 2h) First, write a formula as a linear combination of the function values as where a, b, and c are unknowns to be determined so that we obtain the highest accuracy possible.

Deriving FD Coefficients

LeVeque (2007, p. 7)

17

Rearranging the Taylor’s expansions of u(x – h), and u(x – 2h), we obtain To obtain an approximation of u'(x), we need Thus,

Example 3 (continued)

LeVeque (2007, p. 7)

18

Use the method of undetermined coefficients to determine the coefficients of a finite difference approximation to u"(x) based on the function values u(x – h), u(x), and u(x + h). What is the order of accuracy of the formula?

Exercise

The method of undetermined coefficients can be extended to compute FD coefficients for using points: Assume that u(x) is at least n+1 times differentiable in the interval containing and all stencil points. Taylor series expansions of u at each point xi about are

Deriving FD Coefficients

LeVeque (2007, p. 10)

“We want to find a linear combination of these values that agrees with as well as possible. So we want where p is as large as possible.” It can be shown that the coefficients ci satisfy for i = 1, …, n. This leads to a linear system which has a unique solution.

Deriving FD Coefficients

LeVeque (2007, p. 10)

% k = derivative order % x = location of approximation % xi = set of FD stencil points function c = fdcoeff(k,x,xi) n = length(xi); A = ones(n); xrow = (xi(:)-x)'; for i=2:n A(i,:)=(xrow.^(i-1))./factorial(i-1); end b = zeros(n,1); b(k+1) = 1; c = A\b; c = c(:);

MATLAB Code

LeVeque (2007, p. 11)

▪ “Another way to derive FD coefficients is to approximate the function u(x) by some polynomial p(x) and then use p'(x) as an approximation to u’(x).” (LeVeque 2007, p. 8) ▪ If we want to use the value of function u(x) at n + 1 points, say, xi, i = 0, 1, …, n, p(x) will be a polynomial

• f degree up to n:

where ai(x) are the Lagrange basis functions. ▪ The degree of polynomial is the order of approximation.

Deriving FD Coefficients

Using two points xi and xi+1, the interpolant is and its derivative is which is the forward FD formula.

Example 4

24

The interpolant can be written in terms of Lagrange basis polynomial Li(x) as where

Lagrange Interpolant

Determine the FD coefficients for u"(x) using 3 points by differentiating the corresponding polynomial interpolant.

Exercise

▪ Recall the centered FD approximation ▪ Let’s assume that the problem is periodic, i.e., and . ▪ For each index i we have one equation. ▪ So we have totally n linear equations.

Differentiation Matrix

Trefethen (2000, p. 1)

▪ The linear system can be written in matrix form as ▪ This matrix is Toeplitz and circulant.

Differentiation Matrix

Image from Trefethen (2000, p. 2)

Differentiation matrix

▪ The linear system for a 4th-order FD approximation is

Differentiation Matrix

Image from Trefethen (2000, p. 3)

29

Example 4 (Trefethen, 2000, p. 3): Using 4th-

• rder FD to approximate the first derivative of

using various values of N.

Order of Accuracy (Revisited)

Image from Trefethen (2000, p. 4)

( )

p

E N DN − = ( )

p

E h Ch = h L N =

Using global interpolation (all points are used) will provide the highest-order approximation.

Spectral Method

Image from Trefethen (2000, p. 5)

Using the spectral method, the approximation will converge to exact solution very rapidly until rounding errors prevent any further improvement (Trefethen, 2000).

Spectral Method (continued)

Image from Trefethen (2000, p. 6)

32

Laplace’s equation: Solutions of Laplace’s equation are analytic , i.e., the Cauchy-Riemann equations are satisfied, and called harmonic functions and sum of them is also a solution, i.e., the superposition principle.

Laplace and Poisson Equations

Reference: http://en.wikipedia.org/wiki/Laplace's_equation

Poisson’s equation

33

Consider 2D incompressible, irrotational flow Show that the scalar velocity potential and the stream function defined by satisfied the Laplace’s equation. Vector calculus identify:

Exercise: Fluid Flow

Reference: http://en.wikipedia.org/wiki/Laplace's_equation

34

Given the Maxwell’s equations Show that the scalar potential defined by satisfies the Poisson’s equation

Exercise: Electrostatics

Reference: http://en.wikipedia.org/wiki/Poisson's_equation

35

Given the Maxwell’s equations Show that the vector potential A defined by (Coulomb gauge) satisfies the Poisson’s equation

Exercise: Magnetostatics

Reference: http://farside.ph.utexas.edu/teaching/em/lectures/node40.html

36

Given the Newton’s law of gravitation and Gauss’s law for gravity Show that the gravity potential defined by satisfies the Poisson’s equation

Exercise: Gravity

Reference: http://en.wikipedia.org/wiki/Gauss's_law_for_gravity

37

Given the 1D boundary value problem Using the centered finite difference we obtain the difference equation

1D Poisson Equation

Reference: LeVeque (2007, p. 15)

Dirichlet BC

38

This leads to a system of m linear equations where

1D Poisson Equation

Reference: LeVeque (2007, p. 16)

39

▪ In this case, matrix A is tridiagonal. The tridiagonal linear system can be efficiently solved by the Thomas’ algorithm which is a special case of Gaussian elimination. ▪ A tridiagonal linear system can be written as

Thomas Algorithm

Image from http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm

40

The solution is obtained by back substitution given by

Thomas Algorithm

41

Local truncation error can be determined by substituting the true solution into the difference equation. Although we do not know the true solution, but if we assume a smooth solution, we can use the Taylor’s expansion to obtain

Local Truncation Error

LeVeque (2007, p. 17)

42

Using the Poisson equation we obtain This leads to the linear system where is the vector true solution is the vector of truncation error

Local Truncation Error

LeVeque (2007, p. 17)

43

Global error defined by satisfies the boundary value problem Since , the global error is roughly

Global Error

LeVeque (2007, p. 18)

44

To determine if the FD scheme is stable, consider the linear system with solution with norm Since so the FD scheme is stable if is bounded, i.e.

Stability

LeVeque (2007, p. 18-19)

Triangular inequality

45

Higher-Order Methods

▪ Higher-order finite-difference methods

• Use longer FD stencils.

▪ Richardson extrapolation

• Use low-order FD solutions on two grid

levels to obtain higher-order solution. ▪ Deferred correction

• Use low-order FD solution to

approximate the local truncation error and then solve for the global error to

• btain more accurate solution.

Recall that Let . We than have If we also use another step size , we also have Solving for yields r-order accurate value

Richardson Extrapolation

Heath (2002, p. 369)

47

Richardson Extrapolation

LeVeque (2007, p. 53-54)

Let’s use second-order FD. is true solution. Coarse grid solution Fine grid solution We then have Show that the extrapolated value is 4th-order accurate

48

Deferred Correction

LeVeque (2007, p. 53-54)

First solve to obtain 2nd-order accurate solution. The global error satisfies where We can use the approximate solution to approximate and obtain an estimated global error E to improve the accuracy of u. In this case, we can use the Poisson equation to obtain

49

Fick’s first law: Mass diffusion Concentration goes from regions with high concentration to regions with low concentration where J is diffusion flux D is diffusion coefficient  is substance concentration

Fick’s Laws of Diffusion

http://en.wikipedia.org/wiki/Fick's_laws_of_diffusion

50

A drag or viscosity force exerting on spherical

• bjects with low Reynolds numbers (e.g. small

in size) moving through a viscous fluid is given by where  is dynamic viscosity r is radius of spherical objects v is flow velocity

Stokes’ Law of Drag Force

http://en.wikipedia.org/wiki/Stokes'_law

51

Using Fick’s laws and Stokes’ law, Einstein was able to unify the physical diffusion of Fourier and stochastic diffusion of Laplace through the theory of Brownian motion (Einstein, 1905). Mean squared displacement and time relation: Stokes-Einstein relation:

Stokes-Einstein Relation

http://en.wikipedia.org/wiki/Einstein_relation_(kinetic_theory)

where kB is Boltzmann’s constant T is temperature

52

A diffusion process is governed by the diffusion equation In mass transport, the process is characterized by the diffusion coefficient. In heat conduction, the process is characterized by the diffusion coefficient and thermal conductivity.

Diffusion Equation

http://en.wikipedia.org/wiki/Diffusion_equation

53

Assuming a steady state and constant diffusion coefficient, the diffusion equation becomes the Poisson’s equation which is an elliptic PDE. We can use 2nd-order FD to solve the problem.

Steady-State Problems

http://en.wikipedia.org/wiki/Diffusion_equation

5-point stencil 9-point stencil

• 4

1

• 20

4 1

54

Both stencils are second-order accurate but the 9-point stencil has an advantage when it is used with the deferred correction method.

Accuracy of 2D Stencils

LeVeque (2007, p. 64-65)

55

Consider the 1D time-dependent diffusion equation We can solve the equation using the explicit FD Truncation error: Exercise: Derive the stability condition

1D Diffusion Equation

LeVeque (2007, p. 182)

56

where D2 is the 5-point stencil. This leads to the linear system which must be solved simultaneously to get the approximate solution u. Truncation error:

Crank-Nicolson Method

LeVeque (2007, p. 182)

▪ Einstein, A., 1905, On the motion required by the molecular kinetic theory of heat, of particles suspended in fluids at rest, Annalen der Physik 17, 549-560. ▪ Heath, M.T., 2002, Scientific Computing: An Introductory Survey, McGrawHill. ▪ LeVeque, R. J., 2007, Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems, SIAM. ▪ Trefethen, L. N. 2000, Spectral Methods in MATLAB, SIAM.

References

References

▪ Heath, M. T., 2002. Scientific Computing: An Introductory Survey, Second Edition, McGraw-Hill. ▪ LeVeque, R. J., 2007, Finite Difference Methods for Ordinary and Partial Differential Equation, SIAM. ▪ Trefethen, L. N., 2000, Spectral Methods in MATLAB, SIAM. ▪ Einstein, A., 1905, On the motion required by the molecular kinetic theory of heat, of particles suspended in fluids at rest, Annalen der Physik 17, 549-560.