Atmospheric and Oceanic Modeling - Series Expansion Methods

Lecture 3 Series Expansion Methods Series expansion methods are the general class that encompass spectral and finite element methods. We approximate functions as a linear combination of prescribed expansion functions - we call these basis functions. For a con­ tinuous function f (x), we write f (x) = N ⎤ aj �j (x) (3.1) j=1 where �j , j = 1, . . . , N, are the basis functions that each satisfy any boundary conditions on f (x). The coefficients aj are the unknowns and form a vector of N numbers. Suppose that we have a partial differential equation of the arbitrary form L (f ) = �(x). (3.2) We define the residual of equation (3.2) as ⎪ r(x) = L (f ) − �(x) = L ⎛ N ⎤ j=1 ⎜ aj �j (x)⎞ − �(x). (3.3) If L is linear and the basis functions, �j (x) are the eigen-functions of L, then the residual can be set to zero for the whole domain and the resulting N algebraic equations can be solved for aj . Generally speaking, L is non-linear so we will describe more general approaches for solving for the coefficients aj . We will consider three strategies: i) minimization of the of l2 -norm of the residual, ⎨ 2 (||r(x)||2 ) = r(x)2 dx, (3.4) ii) the collocation method where we set the residual to zero at a discrete set of positions xk (e.g. on a regular grid xk = k�x), r(xk ) = 0 � k = 1, . . . , N (3.5) and iii) the Galerkin method which requires the residual to be orthogonal to each of the basis functions, ⎨ �k r(x)dx = 0 � k = 1, . . . , N (3.6) The collocation method is used in the pseudo-spectral method while the Galerkin method is used extensively in the finite element method. The spec­ tral method is a special case where the l2 -norm and Galerkin method become equivalent. There are many variants on these methods and we will discuss one, the Petrov-Galerkin method. Here, the residual is made orthogonal to a set of test functions, νk (x), which may be different to the basis set, �k (x): ⎨ νk r(x)dx = 0 � k = 1, . . . , N. (3.7) The Petrov-Galerkin method is more general than the Galerkin method be­ cause if we chose νk (x) = �k (x) we recover equation (3.6). 3.1 The Spectral Method Spectral methods are a special case of the series expansion method; the basis functions form an orthogonal set: ⎨ �i �j dx = 0 � i �= j. The ability to use an orthogonal basis set is largely dictated by the domain geometry and boundary conditions. For example, it is natural to use spherical harmonics in spectral atmospheric models but difficult to model irregular coasts with the same representation. 3.1.1 Spectral method compared to the finite differ­ ence method In lecture 2 we analyzed the numerical dispersion of waves using the finite difference method applied to the linear advection problem �t ν + c�x ν = 0 in which c was constant. Consider using a Fourier series expansion to repre­ sent ν(x, t): ν(x, t) = N ⎤ ak (t)eikx k=−N We’re implicitly assuming a periodic domain −∂ � ∂. Note, in this case, ν(x, t) is real so the coefficients in the series satisfy the property ak = a�−k where a�k is the complex conjugate of ak . Direct substitution of each component of the series, ak (t)eikx , into the governing equation yields �t ak + ickak = 0 � k which is unusually trivial to solve. The corresponding dispersion relation shows us that all waves, including the shortest, propagate with the exact correct phase speed. Note, we did not use any of the methods listed in section 3. Let us apply the Galerkin approximation (equation 3.6) to see what happens. Note that two complex functions, g(x) and h(h), are orthogonal is the integral over the domain, S, of the product of one with the complex conjugate of the other is zero: ⎨ g(x)h� (x) dx = 0 S For this problem, the Galerkin approximation is ⎨ � −� ⎭ e−ijx ⎡�t N ⎤ ak (t)eikx + c�x N ⎤ k=−n k=−n ⎣ ak (t)eikx � dx = 0 � j = −N, . . . , N (3.8) and we need to evaluate the following integral: ⎨ � −� e−ijx eikx dx = ⎟ ⎩ i(k−j)x �� e � � = 0 j �= k ⎠ k−j −� � � ⎧ 2∂ j=k Thus the Galerkin approximation yields 2∂�t ak + 2ic∂kak = 0 which is exactly the same results has obtained by direct substitution. The spectral approximation does not introduce phase speed or amplitude errors, ignoring time-discretization errors. 3.1.2 Spectral Stommel model in 1-D This is a somewhat contrived use of the spectral method but allows us to make a direct comparison with finite difference method used in section 1.5. To re-state the problem, we seek solutions to the differential equation: π�xx ρ + �x ρ = −1 with boundary conditions ρ(0, 1) = 0. We will express the solution, ρ(x), in terms of a sin-series: ρ(x) = N ⎤ aj sin (j∂x) (3.9) j=1 since the functions sin (j∂x) all satisfy the boundary conditions. The residual of the governing equation is r(x) = π N ⎤ aj �xx sin (j∂x) + j=1 N ⎤ aj �x sin (j∂x) + 1 j=1 = −π∂ 2 N ⎤ aj j 2 sin (j∂x) + ∂ j=1 N ⎤ aj j cos (j∂x) + 1 (3.10) j=1 Using the Galerkin method (3.6) we require the residual to be orthogonal to the basis functions: ⎨ 0 1 r(x) sin (k∂x) dx = ⎨ 1 0 ⎟ ⎠ −π∂ ⎧ +∂ N ⎤ j=1 2 N ⎤ aj j 2 sin (k∂x) sin (j∂x) j=1 � ⎦ aj j sin (k∂x) cos (j∂x) + sin (k∂x)⎫ dx = 0 � k = 1, . . . , N (3.11) Once we evaluate the integrals and sums it will become apparent that (3.11) represents a set of N algebraic equations in the unknowns, ak . The first term in (3.11) involves the expression ⎨ N 1⎤ 0 j=1 2 aj j sin (k∂x) sin (j∂x) dx = N ⎤ aj j j=1 2 ⎨ 1 0 sin (k∂x) sin (j∂x) dx. Evaluating the integral inside the sum we get ⎨ 1 0 sin (k∂x) sin (j∂x) dx = ⎟ � ⎠ � ⎧ 1 2 j=k 0 j �= k which reflects the orthogonality of the basis set. Evaluating the sum over j, the first term in (3.11) becomes −π∂ 2 N ⎤ aj j 2 j=1 ⎨ 1 0 −π∂ 2 sin (k∂x) sin (j∂x) dx = ak k 2 . 2 The last term is similarly straight forward: 1 − cos (k∂) 1 − (−1)k −1 [cos (k∂x)]10 = = k∂ k∂ k∂ 0 which is 2/k∂ when k is odd and is zero when k is even. The beta term is more complicated; the inner integral evaluates to ⎨ 1 sin (k∂x) dx = ⎨ 1 0 sin (k∂x) cos (j∂x) dx = k − k(−1)(j+k) . ∂(k 2 − j 2 ) When j + k is even, which includes j = k, the numerator and integral is zero. When j + k is odd the expression becomes 2k/∂(k 2 − j 2 ). Substituting all these results back into 3.11 we obtain N ⎤ jk(1 − (−1)(j+k) ) (−1)k − 1 −π∂ 2 k 2 aj = � k = 1, . . . , N ak + k2 − j 2 2 k∂ j=1 (3.12) which represents N algebraic equations for the unknowns ak . We pose this as a linear algebra problem of the form A a = b where ⎪ ⎜ a1 ⎬ . ⎝ ⎬ .. ⎝ ⎬ ⎝ a=⎬ ⎬ a ⎛ k .. . ⎪ ⎬ ⎬ ⎬ ⎝, b = ⎬ ⎝ ⎬ ⎞ ⎛ 2 � .. . (−1)k −1 k� .. . ⎜ ⎝ ⎝ ⎝ ⎝ ⎝ ⎞ and A whose elements are given by Ajk = π∂ 2 k 2 jk(1 − (−1)(j+k) ) ψ ij + 2 k2 − j 2 where the symbol ψij is the Kronecka function; ψij = 1 � i = j and ψij = 0 � i �= j. This is a “full” matrix meaning it is not sparse; many elements are non-zero. Recall that the matrix problem corresponding to the second order finite difference Stommel model was a tri-diagonal matrix which is a lot easier to invert. The basis functions and solution using N = 6 are plotted in Fig. 3.1 and examples of higher truncations given in Fig. 3.2. The error for a range of truncations is plotted as a function of N in Fig. 3.3. Note that the slope of the error curves seems to steepen downward as N increases - this is a strong motivation for using the property of spectral methods; they have incredible convergence properties. Basis functions φ (x) j 1 φj(x) 0.5 0 −0.5 −1 0 0.2 0.6 0.4 0.8 1 0.8 1 0.8 1 Components of solution aj φj(x) 0.4 ajφj(x) 0.2 0 −0.2 −0.4 0 0.2 0.6 0.4 Solution ψ(x) = Σj ajφj(x) 1 ψ(x) 0.5 0 −0.5 0 0.2 0.4 0.6 X Figure 3.1: Top: The basis functions used in the spectral method solution of the Stommel model (N=6). Middle: the components of the solution due to each basis function, aj �j (x). Bottom: the numerical solution and analytical solution. 1 0.8 0.6 ψ(x) 0.4 0.2 0 −0.2 −0.4 0 1 0.8 0.6 0.4 0.2 X Figure 3.2: Spectral solutions to the Stommel problem using different length Fourier series, N = 6, 12, 20, 40, 60, 100, 200, 400. 0 10 L1 L2 L∞ −1 10 −2 Error 10 −3 10 −4 10 −5 10 −6 10 0 10 1 2 10 10 3 10 N Figure 3.3: Convergence of the spectral method. The error, measured by the l1 , l2 and l� norms, is plotted as a function of N . Note the increasing rate of convergence as N is increased. φj (x) φj (x) 1 1 x/Δx j−1 j j+1 x/Δx j−1 j j+1 Figure 3.4: a) The chapeau function and b) the top hat function 3.2 Finite elements using Chapeau functions As an example, consider the problem of constant advection in one-dimension: �t u + c�x u + φ�xx u = 0 We’ll describe the dependent variable, u, using the chapeau basis functions: �j (x) = � 0 1− |x−xj | �x � |x − xj | > �x � |x − xj | � �x (3.13) which is plotted in Fig. 3.4a. Using the Galerkin approximation (3.6), we have: N ⎨ ⎤ � j=1 −� �i (�t uj �j + cuj �x �j − φuj �x x�j ) dx = 0 � i = 1, . . . , N (3.14) In this expression we see products of basis functions and derivatives of basis functions. We can evaluate the integrals of these expressions: ⎨ � 4 �x �i �i dx = 6 −� ⎨ � 1 �x �i �i±1 dx = 6 −� ⎨ � −� � �i �x �i dx = 0 1 2 −� ⎨ � −2 �i �xx �i dx = �x −� ⎨ � ±1 �i �xx �i±1 dx = �x −� ⎨ �i �x �i±1 dx = ± 3.5 ω=πk ω =2s c ) 2 k k 2 ω4=2skck(1+2/3sk ) 2 2 ω6=2skck(1+2/3sk (1+4/5*sk )) 2 ω =2s c /(1−4/3s ) 3 fe k k k 2.5 ωΔx/c 2 F.E. O(Δx4) F.D. O(Δx6) 1.5 F.D. O(Δx4) 1 F.D. O(Δx2) 0.5 0 0 0.1 0.2 0.3 0.4 0.5 kΔx/π 0.6 0.7 0.8 0.9 1 Figure 3.5: Dispersion relations for the O(�x2 ), O(�x4 ) and O(�x6 ) fi­ nite difference methods and the O(�x4 ) finite elements method (using the chapeau basis functions). where the last two integrals were carried out by parts. Equation 3.14 then becomes 1 c φ (ui+1 − ui−1 ) − (ui+1 − 2ui + ui−1 ) = 0 (�t ui−1 + 4�t ui + �t ui+1 ) + 6 2�x �x2 or, returning the finite difference notation, φ c ψi ui − ψii u = 0 Ax �t u + �x2 �x where the averaging operator, Ax is defined: 1 1 Ax u = u + ψii u = (ui−1 + 4ui + ui+1 ) 6 6 The finite element method gives fourth order spatial accuracy in this problem. For comparison, the second order finite difference approximation is c φ �t u + ψi ui − ψii u = 0 �x �x2 and the fourth order finite difference approximation: i c 1 φ 1 ψi u − ψii u − ψii (u − ψii u) = 0 �t u + 2 �x 12 �x 6 which uses a five point stencil. The dispersion relation for the undamped waves (φ = 0) in these approximations are plotted in Fig. 3.5. Note that although the stencil of the finite elements method is only three points (as for the second order finite difference approximation), and that the formal truncation is O(�x4 ), the dispersion relation is far more accurate than even the sixth order finite difference approximation. Finite element Stommel model Note that we can trivially see what the finite element approximation to the Stommel problem would be by simply dropping the time-derivative in the above problem. The resulting discretization is exactly the same as the secondorder finite difference discretization. 3.3 Note on Finite Fourier series Different texts use different apparent representations of finite Fourier series which are in fact equivalent. For a real valued function, �(x), one series representation may be written �(x) = a0 + N ⎤ ak cos (kx) + k=1 N ⎤ bk sin (kx) k=1 which has 2N + 1 degrees of freedom and where all coefficients are real. This representation has the advantage that it is immediately obvious that all terms are real. A more succinct series representation is �(x) = N ⎤ ck eikx k=−N which also has 2N + 1 coefficients but where the coefficients ck are complex. A complex number has two components, real and imaginary, and one may wonder if there are in fact 4N + 2 degrees of freedom. If �(x) is real, then the coefficients must satisfy a certain constraint as we derive now. I R I Let us write each coefficient ck = cR k + ick where ck and ck are real numbers. Then each term in the series is I ck eikx = (cR k + ick ) [cos (kx) + i sin (kx)] = ⎩ � ⎩ � I R I cR k cos (kx) − ck sin (kx) + i ck sin (kx) + ck cos (kx) . It is then obvious that the sum of contributions from modes k and −k can be written in term of cos (kx) and sin (kx) alone: ck eikx + c−k e−ikx = ⎩ R I I (cR k + c−k ) cos (kx) + (c−k − ck ) sin (kx) ⎩ � � R I I +i (cR k + c−k ) sin (kx) + (ck + c−k ) cos (kx) . To ensure that the function has no imaginary component we need to ensure R I I cR k − c−k = 0 and ck + c−k = 0 or more simply ck = c�−k . We can also associate each term in the two forms of Fourier series: a0 = c0 (which is real) R R ak = ck + c−k = 2cR k I I bk = c−k − ck = −2cIk .