This preview shows page 1. Sign up to view the full content.
Unformatted text preview: M.Sc. in Computational Science Fundamentals of Atmospheric Modelling
´ Peter Lynch, Met Eireann
Mathematical Computation Laboratory (Opp. Room 30) Dept. of Maths. Physics, UCD, Belﬁeld. January–April, 2004. Lecture 11 Time Integration Methods 2 Introduction
• We have studied various simple solutions of the shallow water equations by making approximations. • In particular, by means of the perturbation method the equations have been linearised, making them amenable to analytical investigation. • However, to obtain solutions in the general case, it is necessary to solve the full nonlinear system. • In numerical weather prediction (NWP) the fully nonlinear primitive equations are solved by numerical means. • In the atmosphere, the nonlinear advection process is a dominant factor. • To get some idea of the methods used, we look at the simple problem of formulating time-integration algorithms for the solution of the simple advection equation. • The development of accurate and eﬃcient numerical schemes is a huge area of ongoing research, and we can do no more than to introduce the subject. ´ • Much of our research eﬀort in Met Eireann over the past ten years has been in this area, and a brief outline of this work will be given. Discretization Methods
There are several distinct approaches to the formulation of computer methods for solving diﬀerential equations. We will conﬁne ourselves to the ﬁnite diﬀerence methods. Other approaches include ﬁnite element method and the spectral method. The central idea of the ﬁnite diﬀerence approach is to approximate the derivatives in the equation by diﬀerences between adjacent points in space or time, and thereby reduce the diﬀerential equation to a diﬀerence equation. • An analytical problem becomes an algebraic one. • A problem with an inﬁnite degree of freedom is replaced by one with a ﬁnite degree of freedom. • A continuous problem goes over to a discrete one.
5 The Finite Diﬀerence Method
We start by looking at the Taylor expansion of f (x): f (x + ∆x) = f (x) + f (x).∆x + 1 f (x)∆x2 + [O(∆x3)] 2
1 f (x − ∆x) = f (x) − f (x).∆x + 2 f (x)∆x2 + [O(∆x3)] (1) (2) The higher order terms, represented by O(∆x3), become less important as ∆x becomes smaller. We can neglect these and use (1) or (2) to get an approximation for the derivative of f (x) as follows: f (x + ∆x ) − f (x ) + O(∆x) = fF + O(∆x) ∆x f ( x ) − f ( x − ∆x ) f (x ) = + O(∆x) = fB + O(∆x) . ∆x These are called the forward and backward diﬀerences. Keeping only leading terms, we incur errors of order O(∆x). f (x ) = 6 We can do better than this: subtracting (2) from (1) yields the following: f ( x + ∆x ) − f ( x − ∆x ) + O(∆x2) = fC + O(∆x2) f (x ) = 2∆x which is seen to be of order O(∆x2), therefore more accurate for small ∆x. Adding (1) and (2) gives the corresponding expression for the second derivative: f (x + ∆x) − 2f (x) + f (x − ∆x) + O(∆x2) f (x ) = ∆x 2 These centered diﬀerences are of accuracy O(∆x2). We can continue taking more and more terms, but obviously there is a trade-oﬀ between accuracy and eﬃciency. Fourth-order accurate schemes are sometimes used in NWP, but second order accuracy is more popular. Exercise:
Consider the function f (x) = A sin(2πx/L). We know that the derivative is (2π/L)A cos(2πx/L). • Show that a forward diﬀerence approximation gives sin(π ∆x/L) fF (x) = (2π/L)A cos[2π (x + ∆x/2)/L] · π ∆x/L whereas the centered diﬀerence yields: sin(2π ∆x/L) fC (x) = (2π/L)A cos[2πx/L] · 2π ∆x/L • Compare these to the true derivative f (x) and investigate their behaviour for small ∆x. • Demonstrate thus that the centered diﬀerence is of higher order accuracy. Grid Resolution and Accuracy
The size of the gridstep ∆x determines the accuracy of the numerical scheme. For the simple sine function the error depended on the ratio (∆x/L). For synoptic scale waves in the atmosphere a typical value of L is 1000 km. To make the ratio equal to 0.1 we need to have a grid size of about 100 km. This is larger than the typical gridsizes used in operational NWP. The higher the resolution, that is, the smaller the grid-size, the heavier the computational burden. There is a trade-oﬀ between resolution and accuracy. Linear Computational Instability
We consider the equation describing the conservation of a quantity Y (x, t) following the motion of a ﬂuid ﬂow in one space dimension: dY ∂Y ∂Y = 0. ≡ +u dt ∂t ∂x If the velocity is taken to be constant, u = c, or if we linearise ¯ about a mean ﬂow u = c, the equation becomes ∂Y ∂Y +c = 0. ∂t ∂x This simple equation we will call the linear advection equation. It is analogous to a factor of the usual wave equation ∂2 ∂2 − c2 2 ∂t2 ∂x Y= ∂ ∂ +c ∂t ∂x ∂ ∂ −c Y = 0, ∂t ∂x and the general solution is Y = Y (x − ct).
9 10 Since the advection equation is linear, we can construct a general solution from Fourier components Y = a exp[ik (x − ct)] ; k = 2π/L . We take the following initial condition for Y : Y (x, 0) = a exp[ikx] Next, we approximate the diﬀerential equation by a ﬁnite diﬀerence equation using centered diﬀerences for the space and time derivatives. Let the variables x and t be represented by the horizontal and vertical axes. Positive time corresponds to the upper half plane. The initial data occur on the x-axis. The continuous variables are replaced by discrete gridpoints at their integral values and the problem is solved on a ﬁnite diﬀerence grid. Space-Time Grid:
n=5 n=4 n=3 n=2 n=1 n=0 Space axis horizontal Time axis vertical +——–+——–+——–+——–+——–+——–+ I I I I I I I I I I I I I I +——–+——–+——–+——–+——–+——–+ I I I I I I I I I I I I I I +——–+——–+——–+——–+——–+——–+ I I I I I I I I I I I I I I +——–+——–+——–+——–+——–+——–+ I I I I I I I I I I I I I I +——–+——–+——–+——–+——–+——–+ I I I I I I I I I I I I I I +——–+——–+——–+——–+——–+——–+ m=-3 m=-2 m=-1 m=0 m=1 m=2 m=3 We denote the value of Y at a grid point by:
n Y (m∆x, n∆t) = Ym . Inter-dependency of Points
n+1 n n-1 + X + m-1 X O X m + X + m+1 = 0. Then the ﬁnite diﬀerence approximation to the diﬀerential equation may be written as follows:
n n Ym+1 − Ym−1 2∆t n n Ym+1 − Ym−1 +c 2∆x Solving for the value at time (n + 1)∆t gives
n n Ym+1 = Ym−1 − c ∆t n n (Ym+1 − Ym−1) ∆x The value at the time (n + 1)∆t is obtained by adding a term to the value at (n − 1)∆t; the method is known as the leapfrog method because of this jump over the time n∆t. Note the ratio c ∆t λ≡ . ∆x The value of this will be found to be crucial.
13 The evaluation of the equation at point O involves values of n the variable at points X. Solving for Ym+1 thus requires
n Ym−1 , n Ym−1 and n Ym+1 . Note how the leapfrog scheme splits the grid into two independent subgrids.
14 Solving for the value at time (n + 1)∆t gives c ∆t n n n n Ym+1 = Ym−1 − (Ym+1 − Ym−1) ∆x Assuming that we know the solution up to time n∆t, the values at time (n + 1)∆t can be calculated, and the solution advanced one timestep in this way. Then the whole procedure can be repeated to advance the solution to (n + 2)∆t, and so on. Under what conditions does the solution of the ﬁnite diﬀerence equation approximate that of the original diﬀerential equation? Intuitively, we would expect that a good approximation would be obtained provided the grid steps ∆x and ∆t are small enough. However, it turns out that this is not enough, and that the value of the ratio λ = c∆t/∆x is critical. This surprising result has important practical implications. The CFL Stability Criterion
Let us assume a solution of the ﬁnite diﬀerence equation in the form n Ym = a An exp[ikm∆x] Substituting this in the equation and dividing by a common factor gives A2 + (2iλ sin k ∆x)A − 1 = 0 where λ = c∆t/∆x. This is a quadratic equation for the amplitude A, with solutions A± = −iλ sin k ∆x ± 1 − λ2 sin2 k ∆x . Question: The roots of a quadratic equation may be either real or complex, depending on the value of the coeﬃcients. We consider in turn the two possible cases. Case I: |λ| ≤ 1
The quantity under the square-root sign is positive, so the modulus of A is given by |A|2 = 1 − λ2 sin2 k ∆x + (λ sin k ∆x)2 = 1 . The modulus of A is seen to be unity. Thus, we may write A = exp(iψ ) . The two values of the phase are ψ1 = − arcsin(λ sin k ∆x) and ψ2 = π − ψ1 . For small λ, we have ψ1 ≈ −λk ∆x = −kc∆t and ψ2 ≈ π + λk ∆x = π + kc∆t The solution of the equation may now be written
n Ym = D exp(iψ1n) + E exp[i(−ψ1 + π )n] exp(ikm∆x) Taking account of the initial conditions, we get the solution
n Ym = (a − E ) exp[ik (m∆x + ψ1n/k )] + (−1)nE exp[ik (m∆x − ψ1n/k )] Physical Mode Computational Mode Exercise: Check the algebra leading to this solution. The ﬁrst term of this solution is called the physical mode and corresponds to the solution of the diﬀerential equation. The second term is called the computational mode. It arises through the use of centered diﬀerences resulting in the approximation of a ﬁrst order diﬀerential equation by a second order diﬀerence equation (with an extra solution). If the ratio λ is small, the ﬁnite diﬀerence solution is given approximately by Y ≈ a exp[ik (m∆x − cn∆t)] which is just the analytical solution. Thus, we would expect good results from the ﬁnite diﬀerence approximation in this case.
17 18 where D and E are arbitrary constants. The centered diﬀerence approximation cannot be used for the ﬁrst timestep because it would involve values at time t = −∆t which are not known. Instead, we must use another approximation for the ﬁrst step, typically an uncentered forward time step. Thereafter, the leapfrog scheme can be used. In essence, this amounts to speciﬁng another “initial contition”, the computational initial condition, at t = ∆t. The value of this determines the amplitude of the computational mode. It should be chosen to minimize this. The above solution implies
0 Ym = a exp(ikm∆x) 1 Ym = [a exp(iψ1) − 2E cos ψ1] exp(ikm∆x) 1 0 Requiring E = 0, we ﬁnd that Ym = exp(iψ1)Ym. In this simple Case II: |λ| > 1
Recall that the roots of the quadratic are A± = −iλ sin k ∆x ± 1 − λ2 sin2 k ∆x . If |λ| > 1, there will be some wavelengths for which λ2 sin2 k ∆x > 1 . Then the two roots of the quadratic are pure imaginary A = i −λ sin k ∆x ± λ2 sin2 k ∆x − 1 and therefore either |A+| > 1 or |A−| > 1, i.e., the modulus of one of the roots will exceed unity. In that case the amplitude of the corresponding solution will increase without limit as time increases. That is, the amplitude of the solution of the ﬁnite diﬀerence equation will grow without bound for large time. This phenomenon is called computational instability. case, we can eliminate the computational mode. In general, it is much more diﬃcult. In case of computational instability, the solution of the ﬁnite diﬀerence equation cannot possibly resemble the physical solution. The physical solution remains of constant amplitude for all time. The numerical solution grows without limit with time. We thus require that |λ| ≤ 1. The condition for stability is known as the CFL Criterion: Unconditionally Stable Schemes.
´ A large part of the research eﬀort in Met Eireann recently has been devoted to the development of integration schemes which are free of the CFL constraint. The semi-Lagrangian scheme for advection is based on the idea of approximating the Lagrangian form of the time derivative. It is so formulated that the numerical domain of dependence always includes the physical domain of dependence. This necessary condition for stability is satisﬁed automatically by the scheme. The semi-Lagrangian algorithm has enabled us to integrate the primitive equations using a time step of 15 minutes. This can be compared to a typical timestep of 2.5 minutes for conventional schemes. The consequential saving of computation time means that the operational numerical guidance is available to the forecasters much earlier than would otherwise be the case. We discuss semi-Lagrangian schemes in the next lecture.
22 c ∆t ≤1 ∆x
after Courant, Friedichs and Lewy (1928), who ﬁrst published the result. It implies that, if we reﬁne the space grid, that is, decrease ∆x, we must also shorten the time step ∆t. Thus, halving the grid size in a two dimensional domain results in an eightfold increase in computation time.
View Full Document
This note was uploaded on 01/31/2011 for the course MATH 21A taught by Professor Osserman during the Spring '07 term at UC Davis.
- Spring '07