Numerical Differentiation

Numerical Differentiation Introduction These notes provide a basic introduction to numerical differentiation using finite- difference grids. They consider the interplay between truncation error and roundoff error. Finite-difference grids In a finite-difference grid, a region is subdivided into a set of discrete points. The spacing between the points may be uniform or non-uniform. For example, a grid in the x direction, xmin ≤ x ≤ xmax may be written as follows. First we place a series of N+1 nodes numbered from zero to N in this region. The coordinate of the first node, x0 equals xmin. The final grid node, xN = xmax. The spacing between any two grid nodes, xi and xi-1, has the symbol Δxi. These relations are summarized as equation [1]. x0 = xmin xN = xmax xi – xi-1 = Δxi [1] A non-uniform grid, with different spacing between different nodes, is illustrated below. ●---------●------------●------------------●---~ x0 x1 x2 x3 ~-------●----------●------● xN-2 xN-1 xN For a uniform grid, all values of Δxi are the same. In this case, the uniform grid spacing, in a onedimensional problem is usually given the symbol h. I.e., h = x i – xi-1 for all values of i. In ME 309 we will limit our consideration to one-dimensional finite-difference problems. However, advanced courses consider multiple space dimensions discussed below In two space dimensions a grid is required for both the x and y, directions, which results in the following grid and geometry definitions, assuming that there are M+1 grid nodes in the y direction. x0 = xmin y0 = ymjn xN = xmax yM = ymay xi – xi-1 = Δxi yj – yj-1 = Δyj [2] For a three-dimensional transient problem there would be four independent variables: the three space dimensions, x, y and z, and time. Each of these variables would be defined at discrete points, i.e. x0 = xmin y0 = ymjn z0 = zmkn t0 = tmin xN = xmax yM = ymay zK = zmaz tL = tmay xi – xi-1 = Δxi yj – yj-1 = Δyj zk – zk-1 = Δzk tn – tn-1 = Δtn [3] Any dependent variable such as u(x,y,z,t) in a continuous representation would be defined only at discrete grid points in a finite-difference representation. The following notation is used for a one dimensional problem. Numerical methods Page 2 f k  f ( xk ) [4] This notation can be extended to problems in more than one dimension including transient problems. In the most complex case the notation n uijk  u( xi , y j , z k , t n ) is used to denote the value of the dependent at a particular point in the region, (xi, yj, zk, tn) where the variable is defined. Finite-difference Expressions Derived from Taylor Series The Taylor series provides a simple tool for deriving finite-difference approximations. It also gives an indication of the error caused by the finite difference expression. Recall that the Taylor series for a function of one variable, f(x), expanded about some point x = a, is given by the infinite series, f ( x)  f (a)  df 1 d2 f ( x  a)  dx x  a 2! dx 2 ( x - a) 2  xa 1 d3 f 3! dx3 ( x - a)3  .... [5] xa The “x = a” subscript on the derivatives reinforces the fact that these derivatives are evaluated at the expansion point, x = a. We can write the infinite series using a summation notation as follows:  1 dn f n n  0 n ! dx f ( x)   ( x - a )n [6] x a In the equation above, we use the definitions of 0! = 1! = 1 and the definition of the zeroth derivative as the function itself. I.e., d0f/dx0|x=a = f(a). If the series is truncated after some finite number of terms, say m terms, the omitted terms are called the truncation error. These omitted terms are also an infinite series. This is illustrated below.  1 dn f 1 dn f n ( x a )  ( x - a) n  n n n  0 n! dx x  a n  m 1 n! dx x  a Terms used  Truncation error m f ( x)   [7] In this equation the second sum represents the truncation error, ε m, from truncating the series after m terms. m   1 dn f  n n  m 1 n ! dx ( x - a) n [8] xa The theorem of the mean can be used to show that the infinite-series truncation error can be expressed in terms of the first term in the truncation error, that is m  1 d m1 f (m  1) ! dx m1 ( x - a) m1 x  [9] Numerical methods Page 3 Here the subscript, “x = ξ”, on the derivative indicates that this derivative is no longer evaluated at the known point x = a, but is to be evaluated at x = ξ, an unknown point between x and a. Thus, the price we pay for reducing the infinite series for the truncation error to a single term is that we lose the certainty about the point where the derivative is evaluated. In principle, this would allow us to compute a bound on the error by finding the value of ξ, between x and a, that made the error computed by equation [9] a maximum. In practice, we do not usually know the exact functional form, f(x), let alone its (m+1)th derivative. In using Taylor series to derive the basic finite-difference expressions, we start with uniform onedimensional grid spacing. The difference, Δxi, between any two grid points is the same and is given the symbol, h. This uniform grid can be expressed as follows. Δxi = xi – xi-1 = h or xi = x0 + ih for all i = 0,…,N [10] Various increments in x at any point along the grid can be written as follows: xi+1 – xi-1 = xi+2 – xi = 2h xi-1 – xi = xi – xi+1 = –h xi-1 – xi+1 = xi – xi+2 = –2h [11] Using the increments in x defined above and the notation f i = f(xi) the following Taylor series can be written using expansion about the point x = xi to express the values of f at some specific grid points, xi+1 , xi-1 , xi+2 and xi-2. The conventional Taylor series expression for f(x) in equation [5] can be adapted for use in finite differences by writing an expansion equation about a particular grid point, x = xi, to determine the value of f(x) at another grid point, x i+k. From equation [10], we see that xi+k = xi + kh so that f(xi+k) = f(xi + kh). The difference, x, in the independent variable, x, between the evaluation point, xi + kh, and the expansion point, xi, is equal to kh. Using a = xi as the expansion point and kh as x allows us to rewrite equation [5] as shown below. f ( xi  kh)  f ( xi )  df 1 d2 f kh  dx x  xi 2! dx 2 (kh)2 x  xi  1 d3 f 3! dx3 (kh)3  ..... [12] x  xi The next step is to use the notation that f(xi + kh) = fi+k, and the following notation for the nth derivative, evaluated at x = xi. df fi  dx ' d2 f fi  2 dx '' x  xi ... x  xi dn f fi  n dx n [13] x  xi With these notational changes, the Taylor series in equation [12] can be written as follows. f i  k  f i  f i ' kh  f i '' (kh) 2 (kh) 3  f i '''  ..... 2! 3! [14] Finite-difference expressions for various derivatives can be obtained by writing the Taylor series shown above for different values of k, combining the results, and solving for the derivative. The simplest example of this is to use only the series for k = 1. f i 1  f i  f i ' h  f i '' h2 h3  f i '''  ..... 2 6 We can rearrange this equation to solve for the first derivative, f’i; recall that this is the first derivative at the point x = xi. [15] Numerical methods Page 4 fi '  fi 1  fi h h2 f f  fi ''  fi '''  .....  i 1 i  O(h) h 2 6 h [16] The first term to the right of the equal sign gives us a simple expression for the first derivative; it is simply the difference in the function at two points, f(xi+h) – f(xi), divided by h, which is the difference in x between those two points. The remaining terms in the first form of the equation are an infinite series. That infinite series gives us an equation for the error that we would have if we used the simple finite difference expression to evaluate the first derivative. As noted above, we can replace the infinite series for the truncation error by the leading term in that series. Remember that we pay a price for this replacement; we no longer know the point at which the leading term is to be evaluated. Because of this we often write the truncation error as shown in the second equation. Here we use a capital oh followed by the grid size in parentheses. In general the grid size is raised to some power. (Here we have the first power of the grid size, h = h1.) In general we would have the notation, O(h n). This notation tells us how the truncation error depends on the step size. This is an important concept. If the error is proportional to h, cutting h in half would cut the error in half. If the error is proportional to h2, then cutting the step size in half would reduce the error by ¼. When the truncation error is written with this O(h n) notation, we call n the order of the error. In two calculations, with step sizes h1 and h2, we expect the following relation between the truncation errors, ε1 and ε2 for the calculations. h  2   1  2  h1    n [17] We use the approximation sign (≈) rather than the equality sign in this equation because the error term also includes an unknown factor of some higher order derivative, evaluated at some unknown point in the region. The approximation shown in equation [17] would be an equality if this other factor were the same for both step sizes. Another important idea about the order of the error is that an nth order finite-difference expression will give an exact value for the derivative of an nth order polynomial. Because a Taylor series is a polynomial series, it can represent a polynomial exactly if a sufficient number of terms are used. This is illustrated further below. The expression for the first derivative that we derived in equation [16] is said to have a first order error. We can obtain a similar finite difference approximation by writing the general series in equation [14] for k = -1. This gives the following result. f i 1 3 h2 ''' h  fi  fi h  fi  fi  ..... 2 6 ' '' [18] We can rearrange this equation to solve for the first derivative, f’i; recall that this is the first derivative at the point x = xi. 2 f i  f i 1 f  f i 1 '' h ''' h fi   fi  fi  .....  i  O ( h) h 2 6 h ' Here again, as in equation [16], we have a simple finite-difference expression for the first derivative that has a first-order error. The expression in equation [16] is called a forward difference. It gives an approximation to the derivative at point i in terms of values at that point [19] Numerical methods Page 5 and points forward (in the +x direc-tion) of that point. The expression in equation [19] is called a backwards difference for similar reasons. An expression for the first derivative that has a second-order error can be found by subtracting equation [18] from equation [15]. When this is done, terms with even powers of h cancel giving. f i 1  f i 1  2 f i ' h  2 f i ''' h3 h5  2 f i '''''  ..... 6 120 [20] Solving this equation for the first derivative gives the following result. fi '  fi 1  fi 1 h2 h4 f f  f i '''  f i '''''  .....  i 1 i 1  O(h 2 ) 2h 6 120 2h [21] The finite-difference expression for the first derivative in equation [21] is called a central difference. The point at which the derivative is evaluated, xi, is central to the two points (xi+1 and xi-1) at which the function is evaluated. The central difference expression provides a higher order (more accurate) expression for the first derivative as compared to the forward or backward derivatives. There is only a small amount of extra work (a division by 2) in getting this more accurate result. Because of their higher accuracy, central differences are usually preferred in finite difference expressions. Central difference expressions are not possible at the start of end of a boundary. It is possible to get higher order finite difference expressions for such points by using more complex expressions. For example, at the start of a region, x = x0, we can write the Taylor series in equation [14] for the first two points in from the boundary, x1 and x2, expanding around the boundary point, x0. f1  f 0  f 0' h  f 0'' h2 h3  f 0'''  ..... 2 6 3 ( 2 h) 2 ''' ( 2h) f 2  f 0  f ( 2h )  f  f0  ..... 2 6 ' 0 '' 0 [22] [23] These equations can be combined to eliminate the h2 terms. To start, we multiply equation [22] by 4 and subtract it from equation [23].     ( 2h) 2 ( 2h) 3 ( h) 2 ( h) 3 f 2  4 f1   f 0  f 0' (2h)  f 0''  f 0'''  ...  4 f 0  f 0' (h)  f 0''  f 0'''  .. 2 6 2 6     This equation can be simplified as follows f 2  4 f1  3 f 0   f 0' (2h)  4 f 0''' ( h) 3  .... 6 [24] When this equation is solved for the first derivative at the start of the region a second order accurate expression is obtained. f 0'   f 2  4 f1  3 f 0 h2  f 2  4 f1  3 f 0  f 0'''  ....   O( h 2 ) 2h 3 2h [25] Numerical methods Page 6 A similar equation can be found at the end of the region, x = x N, by obtaining the Taylor series expansions about the point x = xN, for the values of f(x) at x = xN-1 and x = xN-2. This derivation parallels the derivation used to obtain equation [25]. The result is shown below. f N'  f N  2  4 f N 1  3 f N h2 f  4 f N 1  3 f N  f N'''  ....  N  2  O( h 2 ) 2h 3 2h [26] Equations [25] and [26] give second-order accurate expressions for the first derivative. The expression in equation [25] is a forward difference; the one in equation [26] is a backwards difference. The evaluation of three expressions for the first derivative is shown in Table 2-1. These are (1) the second-order, central-difference expression from equation [21], (2) the first-order, forwarddifference from equation [16], and (3) the second-order, forward-difference from equation [25]. The first derivative is evaluated for f(x) = e x. For this function, the first derivative, df/dx = ex. Since we know the exact value of the first derivative, we can calculate the error in the finite difference results. In Table 1, the results are computed for three different step sizes: h = 0.4, h = 0.2 and h = 0.1. The table also shows the ratio of the error as the step size is changed. The next-to-last column shows the ratio of the error for h = 0.4 to the error for h = 0.2. The final column shows the ratio of the error for h = 0.2 to the error for h = 0.1. Table 1 Tests of Finite-Difference Formulae to Compute the First Derivative – f(x) = exp(x) x x x x x x x 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.8221 2.0138 2.2255 2.4596 2.7183 3.0042 3.3201 3.6693 4.0552 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.8221 2.0138 2.2255 2.4596 2.7183 3.0042 3.3201 3.6693 4.0552 Results using second-order central differences 1.8221 2.0138 2.0171 0.0034 2.2255 2.2404 0.0149 2.2293 0.0037 2.4596 2.4760 0.0164 2.4637 0.0041 2.7183 2.7914 0.0731 2.7364 0.0182 2.7228 0.0045 3.0042 3.0242 0.0201 3.0092 0.0050 3.3201 3.3423 0.0222 3.3257 0.0055 3.6693 3.6754 0.0061 4.0552 Results using first-order forward differences 1.8221 2.2404 0.4183 2.0171 0.1950 1.9163 0.0942 2.0138 2.4760 0.4623 2.2293 0.2155 2.1179 0.1041 2.2255 2.7364 0.5109 2.4637 0.2382 2.3406 0.1151 2.4596 3.0242 0.5646 2.7228 0.2632 2.5868 0.1272 2.7183 3.3423 0.6240 3.0092 0.2909 2.8588 0.1406 3.0042 3.3257 0.3215 3.1595 0.1553 3.3201 3.6754 0.3553 3.4918 0.1717 3.6693 3.8590 0.1897 4.0552 4.02 2.15 2.15 2.15 2.15 2.15 4.01 4.01 4.01 4.01 4.01 2.07 2.07 2.07 2.07 2.07 2.07 2.07 Numerical methods Page 7 Table 1 Tests of Finite-Difference Formulae to Compute the First Derivative – f(x) = exp(x) h = .4 h = .2 h = .1 Error Ratios Exact x f(x) (h=.4)/ (h=.2)/ f'(x) f’(x) Error f’(x) Error f’(x) Error (h=.2) (h=.1) Results using second-order forward differences 0.6 1.8221 1.8221 1.6895 0.1327 1.7938 0.0283 1.8156 0.0066 4.69 4.32 0.7 2.0138 2.0138 1.9825 0.0313 2.0065 0.0072 4.32 0.8 2.2255 2.2255 2.1910 0.0346 2.2175 0.0080 4.32 0.9 2.4596 2.4596 2.4214 0.0382 2.4508 0.0088 4.32 1.0 2.7183 2.7183 2.6761 0.0422 2.7085 0.0098 4.32 1.1 3.0042 3.0042 2.9934 0.0108 1.2 3.3201 3.3201 3.3082 0.0119 1.3 3.6693 3.6693 1.4 4.0552 4.0552 For the second-order formulae, the error ratios in the last two columns of Table 2-1 are about 4, showing that the second-order error increases by a factor of 4 as the step size is doubled. For the first order expression, these ratios are about 2. This shows that the error increases by the same factor as the step size for the first order expressions. The expected values of the error ratios are only obtained in the limit of very small step sizes. We see that the values in the last column of this table (where the actual values of h are smaller than they are in the next-to-last column) are closer to the ideal error ratio. Truncation errors are not the only kind of error that we encounter in finite difference expressions. As the step sizes get very small the terms in the numerator of the finite difference expressions become very close to each other. We lose significant figure when we do the subtraction. For example, consider the previous problem of finding the numerical derivative of f(x) = e x. Pick x = 1 as the point where we want to evaluate the derivative. With h = 0.1 we have the following data for calculating the derivative by the central-difference formula in equation [21]. f i '  f ' ( x)  f i 1  f i 1 f ( x  h)  f ( x  h) 3.004166  2.722815    2.722815 2h 2h 2(0.1) Since the first derivative of ex is ex, the correct value of the derivative at x = 1 is e1 = 2.718282; so the error in this value of the first derivative is 4.5x10-3. For h = 0.0001, the numerical value of the first derivative is found as follows. f ' ( x)  f ( x  h)  f ( x  h) 2.7185536702  2.7180100139   2.718281832990 2h 2(0.0001) Here, the error is 4.5x10-9. This looks like our second-order error. We cut the step size by a factor of 1,000 and our error decreased by a factor of 1,000,000, as we would expect for a second order error. We are starting to see potential problems in the subtraction of the two numbers in the numerator. Because the first four digits are the same, we have lost four significant figures in doing this subtraction. What happens if we decrease h by a factor of 1,000 again? Here is the result for h = 10-7. f ' ( x)  f ( x  h)  f ( x  h) 2.71828210028724  2.71828155660388   2.7182851763 2h 2(0.0000001) Numerical methods Page 8 Our truncation analysis leads us to expect another factor of one million in the error reduction as we decrease the step size by 1,000. This should give us an error of 4.5x10 -15. However, we find that the actual error is 5.9x10-9. We see the reason for this in the numerator of the finite difference expression. As the difference between f(x+h) and f(x-h) shrinks, we are taking the difference of nearly equal numbers. This kind of error is called roundoff error because it results from the necessity of a computer to round off real numbers to some finite size. (These calculations were done with an excel spreadsheet which has about 15 significant figures.) Figure 2-1 shows the effect of step size on error for a large range of step sizes. For the large step sizes to the right of Figure 2-1, the plot of error versus step size appears to be a straight line on this log-log plot. This is consistent with equation [17]. If we take logs of both sides of that equation and solve for n, we get the following result.   log 2    1   log(  2 )  log(  1 ) n  h  log( h2 )  log( h1 ) log 2   h1  [27] Equation [27] shows that the order of the error is just the slope of a log(error) versus log(h) plot. If we take the slope of the straight-line region on the right of Figure 2-1, we get a value of approximately two for the slope, confirming the second order error for the central difference expression that we are using here. However, we also see that as the step size reaches about 10-5, the error starts to level off and then increase. At very small step sizes the numerator of the finite-difference expression becomes zero on a computer and the error is just the exact value of the derivative. Figure 1. Effect of Step Size on Error 1.E+01 1.E+00 1.E-01 1.E-02 Error 1.E-03 1.E-04 1.E-05 1.E-06 1.E-07 1.E-08 1.E-09 1.E-10 1.E-11 1.E-18 1.E-16 1.E-14 1.E-12 1.E-10 1.E-08 Step Size 1.E-06 1.E-04 1.E-02 1.E+00 Numerical methods Page 9 Final Observations on Finite-Difference Expressions from Taylor Series The notes above have focused on the general approach to the derivation of finite-difference expressions using Taylor series. Such derivations lead to an expression for the truncation error. That error is due to omitting the higher order terms in the Taylor series. We have characterized that truncation error by the power or order of the step size in the first term that is truncated. The truncation error is an important factor in the accuracy of the results. However, we also saw that very small step sizes lead to roundoff errors that can be even larger than truncation errors. The use of Taylor series to derive finite difference expressions can be extended to higher order derivatives and expressions that are more complex, but have a higher order truncation error. The central-difference expression for the second derivative can be found by adding equations [15] and [18]. f i 1  f i 1  2 f i  2 f i '' h2 h4  2 f i ''''  ..... 2 24 [28] We can solve this equation to obtain a finite-difference expression for the second derivative. f i ''  2 f i 1  f i 1  2 f i f  f i1  2 f i '''' h  f  .....  i 1  O( h 2 ) i 2 h 12 h2 [29] Although we have been deriving expressions here for ordinary derivatives, we will apply the same expressions to partial derivatives. For example, the expression in equation [29] for the second derivative could represent d2f/dx2 or ∂2f/∂x2. The Taylor series we have been using here have considered x as the independent variable. However, these expressions can be applied to any coordinate direction or time. Although we have used Taylor series to derive the finite-difference expressions, they could also be derived from interpolating polynomials. In this approach, one uses numerical methods for developing polynomial approximations to functions, then takes the derivatives of the approximating polynomials to approximate the derivatives of the functions. A finite-difference expression with an nth order error that gives the value of any quantity should be able to represent the given quantity exactly for an nth order polynomial.* The expressions that we have considered are for constant step size. It is also possible to write the Taylor series for variable step size and derive finite difference expressions with variable step sizes. Such expressions have lower-order truncation error terms for the same amount of work in computing the finite difference expression. Although accuracy tells us that we should normally prefer central-difference expressions for derivatives, some special applications use one-sided differences. The main example of this is * If a second order polynomial is written as y = a + bx + cx2; its first derivative at a point x = x0 is given by the following equation: [dy/dx]x=x0 = b + 2cx0. If we use the second-order central-difference expression in equation [21] to evaluate the first derivative, we get the same result as shown below: dy y ( x0  h)  y ( x0  h) a  b( x0  h)  c( x0  h) 2  [a  b( x0  h)  c( x0  h) 2 ]   dx x  x0 2h 2h 2 2 2 2 2bh  c( x0  2 x0 h  h )  c( x0  2 x0 h  h ) 2bh  4cx0 h    b  2cx0 2h 2h Numerical methods Page 10 such as convection terms in computational fluid dynamics where one-sided differences, called upwind differences, are used. In solving differential equations by finite-difference methods, the differential equation is replaced by its finite difference equivalent at each node. This gives a set of simultaneous algebraic equations that are solved for the values of the dependent variable at each grid point. Finite difference expressions can be derived from Taylor series. This approach leads to an expression for the truncation error that provides us with knowledge of how this error depends on the step size. This is called the order of the error. In finite-difference approaches, we need to be concerned about both truncation errors and roundoff errors. Roundoff errors were more of a concern in earlier computer applications where limitations on available computer time and memory restricted the size of real words, for practical applications, to 32 bits. This corresponds to the single precision type in Fortran or the Single type in VBA. With modern computers, it is possible to do routine calculations using 64-bit real words. This corresponds to the double precision type in Fortran* or the double type in VBA. The 32-bit real word allows about 7 significant figures; the 64-bit real word allows almost 16 significant figures. * Also known as real(8) or real(KIND=8) in Fortran 90 and later versions; single precision is typed as real, real(4) or real(KIND=4) in these versions of Fortran.