Finite Difference Method (FDM) Chaiwoot Boonyasiriwat August 27, 2020
Introduction to FDM ▪ 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. 2
Finite Difference Approximation ▪ 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.
Commonly Used Formulas ▪ Forward finite-difference scheme ▪ Backward finite-difference scheme ▪ Centered finite-difference scheme
Example 1 D 0 u ( x ) is obviously more accurate than D - u ( x ) and D + u ( x ) LeVeque (2007, p. 4)
Asymmetric Approximation Another finite-difference scheme
Order of Accuracy First-order accurate approximation First-order accurate approximation Second-order accurate approximation Third-order accurate approximation
Example 2 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: LeVeque (2007, p. 4-5)
Example 2 (continued) According to the results, we can see that In general, the approximation error can be written as where p is the order of approximation . LeVeque (2007, p. 5)
Example 2 (continued) Absolute errors versus h on a log-log scale |E| h LeVeque (2007, p. 6)
Truncation Error Taylor’s series expansion is the standard approach for analyzing the truncation error of a finite-difference approximation.
Truncation Error (continued) Rearranging the first expansion yields Truncation error 12
Big O Notation 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 x 0 such that If a truncation error is then 13 https://en.wikipedia.org/wiki/Big_O_notation
Index Notation and Convention To make finite difference expressions concise and simple to handle, we usually use the following index notations. 14
FD Stencil ▪ FD stencil is a graphical representation of a finite difference approximation. ▪ Stencil makes a complicated FD expression looks simpler. -1/12 4/3 -5/2 4/3 -1/12 i - 2 i - 1 i i + 1 i + 2 15
Deriving FD Coefficients ▪ 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 – 2 h ) 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. 16 LeVeque (2007, p. 7)
Example 3 (continued) Rearranging the Taylor’s expansions of u ( x – h ), and u ( x – 2 h ), we obtain To obtain an approximation of u '( x ), we need Thus, 17 LeVeque (2007, p. 7)
Exercise 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? 18
Deriving FD Coefficients 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 x i about are LeVeque (2007, p. 10)
Deriving FD Coefficients “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 c i satisfy for i = 1, …, n . This leads to a linear system which has a unique solution. LeVeque (2007, p. 10)
MATLAB Code % 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(:); LeVeque (2007, p. 11)
Deriving FD Coefficients ▪ “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, x i , i = 0, 1, …, n , p ( x ) will be a polynomial of degree up to n : where a i ( x ) are the Lagrange basis functions. ▪ The degree of polynomial is the order of approximation.
Example 4 Using two points x i and x i +1 , the interpolant is and its derivative is which is the forward FD formula.
Lagrange Interpolant The interpolant can be written in terms of Lagrange basis polynomial L i ( x ) as where 24
Exercise Determine the FD coefficients for u "( x ) using 3 points by differentiating the corresponding polynomial interpolant.
Differentiation Matrix ▪ 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. Trefethen (2000, p. 1)
Differentiation Matrix ▪ The linear system can be written in matrix form as Differentiation matrix ▪ This matrix is Toeplitz and circulant . Image from Trefethen (2000, p. 2)
Differentiation Matrix ▪ The linear system for a 4 th -order FD approximation is Image from Trefethen (2000, p. 3)
Order of Accuracy (Revisited) Example 4 (Trefethen, 2000, p. 3): Using 4 th - order FD to approximate the first derivative of using various values of N. = p E h Ch ( ) = h L N DN − = p E N ( ) 29 Image from Trefethen (2000, p. 4)
Spectral Method Using global interpolation (all points are used) will provide the highest-order approximation. Image from Trefethen (2000, p. 5)
Spectral Method (continued) Using the spectral method, the approximation will converge to exact solution very rapidly until rounding errors prevent any further improvement (Trefethen, 2000). Image from Trefethen (2000, p. 6)
Laplace and Poisson Equations 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 . Poisson’s equation 32 Reference: http://en.wikipedia.org/wiki/Laplace's_equation
Exercise: Fluid Flow 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: 33 Reference: http://en.wikipedia.org/wiki/Laplace's_equation
Exercise: Electrostatics Given the Maxwell’s equations Show that the scalar potential defined by satisfies the Poisson’s equation 34 Reference: http://en.wikipedia.org/wiki/Poisson's_equation
Exercise: Magnetostatics Given the Maxwell’s equations Show that the vector potential A defined by (Coulomb gauge) satisfies the Poisson’s equation 35 Reference: http://farside.ph.utexas.edu/teaching/em/lectures/node40.html
Exercise: Gravity 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 36 Reference: http://en.wikipedia.org/wiki/Gauss's_law_for_gravity
1D Poisson Equation Given the 1D boundary value problem Dirichlet BC Using the centered finite difference we obtain the difference equation 37 Reference: LeVeque (2007, p. 15)
1D Poisson Equation This leads to a system of m linear equations where 38 Reference: LeVeque (2007, p. 16)
Thomas Algorithm ▪ 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 39 Image from http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
Thomas Algorithm The solution is obtained by back substitution given by 40
Local Truncation Error 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 41 LeVeque (2007, p. 17)
Local Truncation Error Using the Poisson equation we obtain This leads to the linear system where is the vector true solution is the vector of truncation error 42 LeVeque (2007, p. 17)
Recommend
More recommend
