Atmospheric and Oceanic Modeling - The Finite Difference Method

Lecture 1 Taylor series and finite differences To numerically solve continuous differential equations we must first define how to represent a continuous function by a finite set of numbers, fj with j = 1, 2, ..., and then translate the continuous differential equations into a finite set of algebraic equations that determine the values, fj . There are two basic strategies: (i) the grid-point methods (finite difference and finite volume) and (ii) series expansion methods (spectral, finite element, ...). Of grid-point methods, the finite difference method is the most widely known in ocean modeling and weather forecasting but in recent years, the finite volume method has become more popular. Spectral methods are widely used in meteorology because they are both accurate and efficient but they are very difficult to use in irregular domains and hence are rarely considered in oceanography. 1.1 Representation of functions In Fig. 1.1, a continuous and differentiable function, f , of the independent variable x, is approximated by the values fi = f (xi ) at the points, xi . The set of points, xi , make up the numerical mesh or “grid”. The spacing between the grid points need not be constant; if the spacing is constant then the grid is regular. Technically, no assumption is made about the value of the approximate solution between the grid points. However, in the space between discrete 1 2 f(x) fi+1 fi fi−1 x xi−1 xi xi+1 Figure 1.1: An approximation to a smooth and continuous function, f (x), can be represented by the truncated set of values {fi−1 , fi , . . .} evaluated at the grid points {xi−1 , xi , . . .} where fi � f (xi ). grid points, for example xi−1 < x� < xi , the function can be interpolated assuming a particular representation of the function based on the known values fi−1 , fi , . . .. For example, Fig. 1.2 shows a straight line drawn between values fi−1 and fi which is given by x� − xi−1 f˜(x� ) = (1−�(x� ))fi−1 +�(x� )fi where �(x) � xi − xi−1 ≤xi−1 � x� � xi (1.1) � ˜ so that f (x ) is the approximation to f (x). This is “linear interpolation”. Higher order interpolation and other representations (such as quadratic, spline or high order polynomial) are possible. 1.2 Taylor series and truncation error Taylor series can be used to estimate a truncation error in an approximation of a function and also used to establish approximations to derivatives. The Taylor series expansion at a point x is f (x ± β) = f (x) + � � j=1 βj (±1)j � j f j! �xj 3 f(x) Δx Δx1 Δx2 fi f(x’) fi−1 x xi−1 x’ xi Figure 1.2: Linear interpolation (dashed line) between values fi of a function f (x) (thin line) on a discrete grid xi . = f (x) ± βf � (x) + β3 β4 β2 �� f (x) ± f ��� (x) + f ���� (x) + . . . 2! 3! 4! where β is some small displacement from x. To estimate the error incurred evaluating the function f (x) using linear interpolation (equation 1.1) we can replace the values fi−1 and fi with Taylor series expansions about the point of evaluation, x: �x21 �� f (x) − . . . 2! �x22 �� fi = f (x + �x2 ) = f (x) + �x2 f � (x) + f (x) + . . . . 2! fi−1 = f (x − �x1 ) = f (x) − �x1 f � (x) + where �x1 � x − xi−1 and �x2 � xi − x. Substituting into (1.1) gives �x2 f (x) − �x1 f � (x) + �x1 + �x2 ⎠ �x1 + f (x) + �x2 f � (x) + �x1 + �x2 1 = f (x) + �x1 �x2 f �� (x) + . . . . 2 f˜(x) = ⎠ 1 2 �� �x f (x) − . . . 2 1 � 1 2 �� �x f (x) + . . . 2 2 � 4 This allows us to estimate the error f˜(x) − f (x) as 1 1 f˜(x) − f (x) � �x1 (�x − �x1 )f �� (x) = �x2 �(1 − �)f �� (x) (1.2) 2 2 dropping higher order terms. Here �x = xi − xi−1 = �x1 + �x2 is the grid spacing and �x1 �x1 �= = �x �x1 + �x2 is the non-dimensional position within the space interval (0 � � � 1). Note some properties about the estimate of error for linear interpolation: • the error scales with the square of the grid spacing (�x2 ), • the error is reduced if the interpolation point is near the ends of the interval (� 2 ≡ 0, 1). i.e. near a data value, • the error is largest in the middle of the interval (� ≡ 21 ), equidistant from the data values. The order of accuracy of an approximation refers to the scaling with mesh parameters (normally the grid resolution). The above analysis tells us that linear interpolation is second order accurate or that errors are of order O(�x2 ). This means that doubling the resolution (i.e. halving the grid spacing) reduces the error by a factor of four. 1.3 Constructing difference operators using Taylor series Finite differencing is the method of approximating partial derivatives with differences between function values (on a grid). The Taylor series for fi−1 (x), fi (x), ... can be re-arranged to yield expressions for the derivatives, f � (x), f �� (x), ... yielding formula for approximating derivatives with discrete differences and truncation error. As an example, consider how to estimate the first derivative, f � , at xi in terms of fi−1 and fi . The Taylor series for fi−1 and fi expanded about xi are: �x2 �� fi−1 = f (xi ) − �xf � (xi ) + f (xi ) − . . . 2! fi = f (xi ) 5 where �x = xi − xi−1 is the grid spacing. Solving for f � (xi ) gives: f � (xi ) � 1 �x �� (fi − fi−1 ) + f (xi ) �x 2 (1.3) where we have dropped higher order terms. The difference formula (1.3) is known as a side difference and is first order accurate, i.e. the truncation error is of order O(�x). This is the simplest and least accurate of the finite difference approximations but will appear in latter discussions about the discretization of advection terms. Now, consider the same data points but this time evaluate the derivative at the mid-point, x¯ = 12 (xi−1 + xi ). The Taylor series are: �x � f (¯ x) 2 ⎞ ⎟ 1 �x 2 �� + f (x̄) − 2! 2 �x � = f (¯ x) + f (¯ x) 2 ⎞ ⎟ 1 �x 2 �� + f (x̄) + 2! 2 x) − fi−1 = f (¯ fi 1 �x 3! 2 ⎞ ⎟3 f ��� (¯ x) + 1 �x 4! 2 ⎞ ⎟4 f ���� (¯ x) − . . . 1 �x 3! 2 ⎟3 f ��� (¯ x) + 1 �x 4! 2 ⎟4 f ���� (¯ x) + . . . ⎞ ⎞ Solving for f � (x̄) gives: f � (x̄) � 1 �x2 ��� (fi − fi−1 ) − x) f (¯ �x 4 · 3! (1.4) This is a centered difference and is second order accurate because the trun­ cation error is of order O(�x2 ). Because the evaluation point is at the mid-point between grid nodes, the derivative is said to be staggered in space. Comparing the formula (1.3) and (1.4), we see that the difference operator is the same in each but the point at which the derivative is being approxi­ mated matters. Although the centered difference is formally of higher order the truncation terms can not be directly compared since f �� (xi ) and f ��� (x̄) are different and not known. As a rule, centered evaluation of difference formula are more accurate than off-center or side differences. Further, centered differences are generally of even order and side differences typically of odd order, though it is possible to construct even order side differences by placing the all points of the stencil to one side of the evaluation point. 6 Finally, consider a centered difference evaluated at a grid point (rather than a mid-point). Limiting ourselves to a uniform grid (i.e. xi+1 − xi = xi − xi−1 = �x) we can construct the operator by applying the formula (1.4) over a 2�x interval, yielding: f � (xi ) � 1 �x2 ��� (fi+1 − fi−1 ) − f (xi ) 2�x 3! (1.5) This is also of second order accuracy but here the coefficient in the truncation error is 4 times larger than for the staggered center difference. Therefore, although we can not directly compare f ��� (xi ) with f ��� (x̄) because they are at different locations, we can infer that the staggered difference is more accurate than the 2�x difference when the variations in the function are well resolved. 1.3.1 Second derivative To find a formula for the second derivative, f �� , we must use a stencil of at least 3 points. Here, we’ll consider a uniform grid with spacing �x. The Taylor series expanded about xi are: fi−1 = f (xi ) − �xf � (xi ) + �x2 �� �x4 ���� �x3 ��� f (xi ) + f (xi ) − f (xi ) − . . . 2! 3! 4! fi = f (xi ) fi+1 = f (xi ) + �xf � (xi ) + �x3 ��� �x2 �� �x4 ���� f (xi ) + f (xi ) + f (xi ) + . . . 2! 3! 4! Solving for f �� (xi ) gives: fi+1 − 2fi + fi−1 2�x2 ���� f (xi ) � − f (xi ) 4! �x2 �� (1.6) dropping higher order terms. 1.3.2 Generating higher order differences The process of finding difference formula and truncation terms has involves several expansions of the Taylor series about the same point and then solving the resulting simultaneous equations for the unknown derivatives. This can be cast as a linear algebra problem and where we use Gaussian elimination to solve for the derivatives. 7 The three Taylor series used in section 1.3.1 can be written as a matrix problem: � 1 −�x � 0 � 1 1 �x �x2 2 − �x 3! 0 0 �x2 2 3 � �x4 4! 0 �x3 3! �x4 4! ⎛ f (x ) i ⎛� � f � (xi ) � ⎜� ⎝ � f �� (xi ) � ��� � f (xi ) f ���� (xi ) � ⎛ ⎜ ⎜ f i−1 ⎜ � ⎜ ⎜ ⎜ = � fi ⎝ ⎜ fi+1 ⎝ where the vertical line indicates a square sub-matrix. If we multiply through by the inverse of this sub-matrix then we get: 1 0 0 � � 0 1 0 0 0 1 � 0 �x2 6 0 � f (xi ) f � (xi ) f �� (xi ) f ��� (xi ) f ���� (xi ) � 0 � ⎜� 0 ⎝� � ⎛ �x2 � � 12 ⎛ � ⎜ ⎜ ⎜ � ⎜ ⎜=� ⎜ ⎝ 0 −1 �x 1 �x2 1 0 −2 �x2 0 1 �x 1 �x2 ⎛� ⎛ fi−1 ⎜� ⎜ ⎝ � fi ⎝ fi+1 This approach is identical to the derivations used in previous sections but in just one operation yields all the difference approximations of derivatives up to the order of the number of points. Using five points, centered at xi , the result is: � 1 0 0 0 0 0 30 � � � =� � � � −�x4 90 0 �x2 4 0 �x2 6 0 0 0 0 0 1 � 0 0 −�x4 � � 0 1 0 0 0 � � � 0 0 1 0 0 � � 0 0 0 1 0 � ⎛ 0 0 1 12�x −1 12�x2 −1 2�x3 1 12�x4 −2 3�x 4 3�x2 1 �x3 −4 �x4 � � � ⎜� ⎜� ⎜� ⎜� ⎜� ⎜� ⎜� ⎝� � � 1 0 −5 2�x2 0 6 �x4 f (xi ) f� f �� f ��� f ���� f ����� f ������ 0 2 3�x 4 3�x2 −1 �x3 −4 �x4 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 −1 12�x −1 12�x2 1 2�x3 1 �x4 ⎛� ⎜� ⎜� ⎜� ⎜� ⎜� ⎜� ⎝� fi−2 fi−1 fi fi+1 fi+2 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ This gives us fourth order accurate formula for f � (x) and f �� (x) and second order accurate formula for f ��� (x) and f ���� (x). 8 1.4 The size of truncation errors It is worth bearing in mind that implicit in deriving the above formula is the assumption that the truncation errors are small. We can formalize this by non-dimensionalizing the problem. Consider the staggered difference given by equation (1.4). If variations in f (x) are of order F and have length scale L then we can non-dimensionalize (1.4): � F �� � F � � F �x� 2 ���� � � f (¯ x )� f − f − f (¯ x ) i−1 L L�x� i 4 · 3!L or � 1 � � �x� 2 ���� � � − f (¯ f − f x ) i−1 �x� i 4 · 3! where �x� = �x/L is a non-dimensional resolution parameter. The actual error, normalized, is: 1 (fi − fi−1 ) −1 βa = �x � f (x) f �� (¯ x� ) � and an estimate of the error can be based on the non-dimensional truncation term: � �⎛ � 1 � � ⎞ ⎟ f − f 1 �x 2 i−1 i �x ⎝−1≡ βe = O � 4 · 3! L f � � (x� ) Example: Using the function f (x) = ex , we can verify that the actual error converges with the right power as �x ≈ 0 and that the estimate of error reflects the actual error. By substituting the continuous solution into the difference approximation we can calculate the actual error, which is: βa = 1 �x/2 (e − e−�x/2 ) − 1 �x while the estimate of the error, based on the the truncation terms is: βe = �x2 24 These errors are plotted in Fig. 1.3 and shows that the estimate is represen­ tative of actual errors for all reasonable �x (i.e. �x < 1). 9 2 10 ε : f(x)=ex/L a ε : f(x)=sin(x/L) a 2 ε =(Δx/L) /24 1 10 e 0 Normalized error (ε) 10 −1 10 −2 10 −3 10 −4 10 −1 10 0 10 Δx/L 1 10 Figure 1.3: Actual and estimated errors for centered finite difference approx­ imation applied to the functions ex/L and sin x/L. For small �x, all curves have a slope of −2. Note, also, that the estimate of error, βe , is reasonable even up to �x ≡ 1. 10 Example: The same analysis can be made for the function f (x) = sin x/L. The actual error, in this case, is: βa = 2L �x sin −1 �x 2L and is also plotted in Fig. 1.3. 1.5 Stommel model in 1-D We can illustrate the process of discretization and analysis with a simple one dimensional version of the Stommel model. We will return to the Stommel model at a later stage and explain where and what it is. The equation to be solved, written in non-dimensional form, is: β�xx γ + �x γ = −1 (1.7) with boundary conditions γ(x = 0) = 0 and γ(x = 1) = 0. Here, β is a nondimensional friction which determines a boundary layer width. An analytic solution to (1.7) exists which is: � x � γ = C e− � − 1 − x with 1 C −1 = e− � − 1. The existence of an analytic solution allows us to compare the estimated ac­ curacy of the discrete equations with the true error of the numerical solution. We will find approximate numerical solutions to (1.7) on a regular grid of N + 1 points, xi = �x(i− 1) where �x = 1/N using second order differences. The second order approximations to the first and second derivatives of a function are given by (1.5) and (1.6). Substituting into (1.7) we obtain the following discrete system of N + 1 simultaneous algebraic equations: β 1 (γi−1 − 2γi + γi+1 ) + (−γi−1 + γi+1 ) = −1 2 2�x �x γi = 0 This can be posed as a linear algebra problem: Aγ=b ≤ i = 2...N ≤ i = 1, N + 1. 11 where A is a tri-diagonal (sparse) matrix given by: Ai,i−1 = Ai,i = Ai,i+1 = A1,1 = AN +1,N +1 = β 1 − 2 �x �x −2β �x2 β 1 + 2 �x �x 1 1 ≤ i = 2...N The tri-diagonal structure of A means it can be solved “directly” using, for example, LU decomposition or cyclic reduction. Any row of the matrix between i = 2 and i = N corresponds to the stencil of the difference equation; ⎠ β 1 − 2 2�x �x −2β �x2 β 1 + 2 �x 2�x � Note that for large β, the signs of the elements of the stencil are [1 − 1 1] and that at a critical value of β = �x/2 the sign of the western most point changes (to the left at i − 1) i.e. the stencil has signs [−1 − 1 1]. When this happens, the equations are no longer elliptic (the hyperbolic term dom­ inates everywhere) and some of the eigenvalues have imaginary components. Physically, this change in behavior occurs when the boundary is on the verge of being resolved. When under resolved, the imaginary components of the eigenvalues lead to oscillations in the solution. Numerical and analytic solutions are plotted in Fig. 1.4. For the well re­ solved case, the two curves are barely indistinguishable, but as the resolution is reduced, significant differences (errors) appear in the boundary layer. In the under-resolved case, oscillations emanate from the boundary layer and reach into the interior. Clearly, If we were more concerned with obtaining a solution that was “robust”, meaning that it’s properties are not critically dependent on numerical pa­ rameters, then we could use first order differencing for the beta term. In this case, the stencil is: ⎠ β �x2 −2β 1 − 2 �x 2�x β 1 + 2 �x �x � 12 Note, the the first order difference is “up-wind” in the sense of westward propagating characteristics. Here, the signs of the stencil elements are inde­ pendent of the numerical parameters and we can then expect the eigenvalues of the problem to remain real also. Numerical solutions are also plotted in Fig. 1.4 and here the solutions all look “physical” (i.e. no non-physical os­ cillations) for all resolutions. However, the convergence properties are poor, formerly of first order and so even when the western boundary layer is well resolved, the solution is visibly inaccurate in the vicinity of the boundary. 1.5.1 Measures of error In this particular example, we can measure the error formally since we have an analytical solution. Williamson, 1992, defined a set of “normalized global errors” for use in evaluating global meteorological models. In this context, they are: l1 (γ) = I(|γ − γT |) I(|γT |) (1.8) 1 l2 (γ) = I(|γ − γT |2 ) 2 1 I(|γT |2 ) 2 max(|γ − γT |) l� (γ) = max(|γT |) (1.9) (1.10) where I(γ) = 01 γ dx and γT is the true solution. The normalized errors for the numerical solutions to the Stommel model are plotted as a function of resolution in Fig. 1.5. The slope confirms that the solutions converge with second and first order accuracy respectively for the centered and upwind schemes. Note that the convergence behaviour is non-trivial at very low resolutions and that in terms of these absolute measures, the first order scheme can appear more accurate than the second order scheme. � 1.6 Common notation and stencils Writing out the algebraic equations using subscripts resulting from differ­ enced model equations can be a tedious and very lengthy process. A no­ tation, used widely since Arakawa and Lamb, 1968, greatly simplifies the process. 13 There are two basic operators, the (staggered) center difference and the average or center linear interpolation. We define the operators as follows: αi f = fi+ 1 − fi− 1 2 2 1 i (f 1 + fi− 1 ) f = 2 2 i+ 2 (1.11) (1.12) These operators satisfy the following rules: 1. The difference and average operators commute: α i αj f = α j αi f αi f f i j = αi j = f (1.13) j j (1.14) i (1.15) 2. Differencing products obeys a discrete analog of the product rule: i αi (f g) = f αi g + g i αi f i αi (f g) = f αi g + gαi f i (1.16) (1.17) 3. Averaging products obeys a sum of squares rule: fg i fg i i 1 i = f g i + (αi f )(αi g) 4 1 i = f g + αi (gαi f ) 4 (1.18) (1.19) These relations and the notation greatly simplify some analysis. For example, to prove that a particular discretization of advection conserves the volume integral of the second moment (< ψu�x ψ >= 0): ψuαi ψ i i 1 i = ψ uαi ψ − αi (u(αi ψ)(αi ψ)) 4 i 2 ψ u − αi ( (αi ψ)2 ) = uαi 2 4 All the centered difference approximations can be constructed from the difference and average operators. For example, the staggered center difference is: fi+ 1 − fi− 1 1 2 2 f� � = αi f �x �x 14 The 2�x difference is: fi+1 − fi−1 1 i = αi f 2�x �x The second order, second derivative is: f� � fi+1 − 2fi + fi−1 1 α i αi f = �x2 �x2 It is sometimes useful to describe the “stencil”, being the pattern of con­ nections in a difference equation or operator. In these notes, we will use square brackets to indicate stencil notation. For example, the second order difference fi+1 − fi−1 i αi f = 2�x has a two point stencil, which using our short-hand is f �� � 1 [−1 0 1] 2�x and the fourth order difference 1 2 2 1 fi−2 − fi−1 + fi+1 − fi+2 f� � 12�x 3�x 3�x 12�x has the five point stencil 1 [1 − 8 0 8 − 1] 12�x We can de-compose this stencil into the difference of two stencils: [1 − 8 0 8 − 1] = [1 − 8 1 0 0] − [0 0 1 − 8 1] and each of these can be de-composed further: [1 − 8 1] = −[0 6 0] + [1 − 2 1]. This tells us that the fourth order finite difference can be written succinctly as: i 1 1 � αi (f − αii f ) . f � �x 6 2 Note that the second term is a discrete approximation of �6x f ��� which is the truncation term of the second order finite difference. In other words, higher order difference formula can be found by substituting a difference approximation for the leading truncation terms in low order formula. 15 n=8 Δx/2ε=3.125 1 Analytical solution n Sub−sampled analytic sol 2 Numerical solution O(Δx ) Numerical solution O(Δx) 0.8 Ψ 0.6 0.4 0.2 0 0 0.2 0.4 n=25 0.6 X Δx/2ε=1 0.8 1 1 Analytical solution n Sub−sampled analytic sol 2 Numerical solution O(Δx ) Numerical solution O(Δx) 0.8 Ψ 0.6 0.4 0.2 0 0 0.2 0.4 n=50 0.6 X Δx/2ε=0.5 0.8 1 1 Analytical solution n Sub−sampled analytic sol Numerical solution O(Δx2) Numerical solution O(Δx) 0.8 Ψ 0.6 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 X Figure 1.4: Solutions to the 1-D Stommel model for three resolutions: underresolved, critical and resolved. Plotted are the analytic solution, the analytic solution sub-sampled onto the model grid and two numerical solutions, one using second order (centered) differences and the other usign first order up­ wind difference of the ∂-term. The oscillations in the second order solution (top panel) result from the appearance of imaginary eigenvalues in the ma­ trix problem. The first order solution is “robust” in that it never exhibits oscillations but it converges more slowly than the second order solution. 16 0 10 Δx1 −1 10 −2 Error 10 Δx2 −3 10 L Centered ∞ L Centered 2 L Centered 1 L Upwind ∞ L Upwind 2 L Upwind 1 −4 10 −5 10 −3 10 −2 10 −1 Δx 10 0 10 Figure 1.5: The convergence of the numerical solutions to the 1-D Stommel model. The errors are the l1 , l2 and l� norms and are plotted for the centered (×) and upwind (∀) differenced models for different resolutions. The slope of the centered difference error curve is 2 which is consistent with the second order accuracy of the discrete equations. Similarly the upwind differenced error curves has a slope of 1 but note the tapering at low resolution.