Properties of numerical methods The following criteria are crucial to the performance of a nume rical algorithm PDF document - DocSlides

Page 1
Properties of numerical methods The following criteria are crucial to the performance of a nume rical algorithm: 1. Consistency The discretization of a PDE should become exac t as the mesh size tends to zero (truncation error should vanish) 2. Stability Numerical errors which are generated during the solution of discretized equations should not be magnified 3. Convergence The numerical solution should approach the ex act solution of the PDE and converge to it as the mesh size tends to zero 4. Conservation Underlying conservation laws should be res pected at the discrete level (artificial sources/sinks are to be avoided) 5. Boundedness Quantities like densities, temperatures, co ncentrations etc. should remain nonnegative and free of spurious wiggles These properties must be verified for each (component of the) n umerical scheme
Page 2
Consistency Relationship: discretized equation differential equation Truncation errors should vanish as the mesh size and time step t end to zero Example. Pure convection equation ∂u ∂t ∂u ∂x = 0 discretized by CDS in space, FE in time: +1 +1 2 [( ( Taylor series expansions: +1 + ∂u ∂t ( ∂t . . . ∂u ∂x ( ∂x ( ∂x . . . Hence, +1 +1 2 ∂u ∂t ∂u ∂x = 0 where ∂t ( ∂x [( ( residual of the difference scheme for the exact nodal values x, m
Page 3
Stability Relationship: numerical solution of discretized equations exact solution of discretized equations Definition 1 Numerical errors (roundoff due to final precision o f computers) should not be allowed to grow unboundedly Definition 2 The numerical solution itself should remain unifo rmly bounded Stability analysis: can only be performed for a very limited ra nge problems Matrix method: Au +1 Bu +1 Cu where is assumed to be a linear operator. In practice = so that +1 Cu for the numerical solution of the discretized equations +1 for the exact solution of the discretized equations +1 Ce for the roundoff error incurred in the solution process
Page 4
Matrix method for stability analysis In the linear case +1 Cu . . . , e +1 Ce . . . i.e., the error evolves in the same way as the solution and is bo unded by || ||≤|| || || || || || ) = max spectral radius of Unstable schemes: if 1 then || || 1 and the errors may grow Example. Convection-diffusion equation ∂u ∂t ∂u ∂x ∂x discretized by CDS in space, FE in time: +1 +1 2 +1 ( or +1 +1 ) + +1 where is the Courant number ( is the diffusion number
Page 5
Matrix method for stability analysis a b c a b c a b c = 1 lo frequencies high frequencies Eigenvectors = cos sin i , , i Eigenvalues a +1 b c divide by i = 1 + 2 (cos 1) + sin = 1 , . . . , N Stability condition = [1 + 2 (cos 1)] + 4 sin pure convection: = 0 ⇒ | | 1 unconditionally unstable :-( pure diffusion: = 0 ⇒ | | 1 if conditionally stable
Page 6
Von Neumann’s stability analysis Objective: to investigate the propagation and amplification of numerical errors Assumptions: linear PDE, constant coefficients, periodic boun dary conditions Continuous error representation x, t ) = ik Fourier series , e ik = cos sin i.e. the error is a superposition of harmonics characterized by their wave number (for wave length ) and amplitide +1 min = 2 max = 2 Discretization m , Here is the phase angle, is the number of waves fitted into the interval L, L ) and determines the highest frequency resolvable on the mesh
Page 7
Von Neumann’s stability analysis Representation of numerical error (trigonometric interpolat ion) i i where Due to linearity, the error satisfies the discretized equati on and so does each harmonic. Hence, it suffices to check stability for i Amplification factor +1 the enhancement of the th harmonic during one time step Stability condition | guarantees that the error component = ( remains bounded Remark. The accuracy of approximation can be assessed by analyzing ph ase errors i.e. the actual speed of harmonics as compared to the exa ct speed
Page 8
Example: pure convection equation in 1D 1. Let ∂u ∂t ∂u ∂x = 0 be discretized by CDS in space, FE in time +1 +1 2 = 0 and substitute iθj The resulting difference equation for the error can be writte n as +1 iθj i +1) i 1) ) = 0 , Divide by iθj +1 i i ) and note that i i = cos sin cos sin = 2 sin Amplification factor +1 = 1 i sin is responsible for stability = 1 + sin the scheme is unconditionally unstable :-(
Page 9
Example: pure convection equation in 1D 2. Let ∂u ∂t ∂u ∂x = 0 be discretized by BDS in space, FE in time +1 = 0 and substitute iθj which yields ( +1 iθj νa iθj i 1) ) = 0 , = 1 νe i = 1 (cos sin ) = 1 sin i sin Re ) = 1 cos θ, Im ) = sin Stability restriction 1 means that must lie within the unit circle in the complex plane. This leads to the CFL condition v > 0 upwind scheme, stable for v < 0 downwind scheme, unconditionally unstable Im Re region of instabilit stabilit region of
Page 10
Example: pure convection equation in 1D The numerical domain of dependence should contain the analyti cal one: if ν > 1, then the data at some grid point may affect the true solution but not the numerical one on the other hand, for ν < 1 some grid points in- fluence the solution although they should not for accuracy reasons it is desirable to have 1; some schemes are exact for = 1 (unit CFL property)















































































































































































+1 +1 dx dt dx dt 3. Let ∂u ∂t ∂u ∂x = 0 be discretized by CDS in space, BE in time +1 +1 +1 +1 2 = 0 and substitute iθj +1 iθj +1 i +1) i 1) ) = 0 1+ i sin It follows that 1+ sin 1 unconditional stability
Page 11
Spectral analysis of numerical errors Consider ∂u ∂t ∂u ∂x ∂x convection-diffusion equation Exact solution: x, t ) = ik vt dt ikx , a ) = ikvt dt a wave with exponentially decaying amplitude traveling at con stant speed Amplification factor ex +1 ikv )( +1) ikv i ex , t, kv arg( ex Amplitude error num ex num numerical damping Phase error arg num arg ex arg num numerical dispersion where arg num is the numerical propagation speed Harmonics travel too fast if 1 (leading phase error) and too slow in the case 1 (lagging phase error)
Page 12
Convergence Relationship: numerical solution of discretized equations exact solution of the differential equation Definition: A numerical scheme is said to be convergent if it pro duces the exact solution of the underlying PDE in the limit 0 and Lax equivalence theorem: stability + consistency = convergence For practical purposes, convergence can be investigated nu merically by comparing the results computed on a series of successively re ned grids The rate of convergence is governed by the leading truncatio n error of the discretization scheme and can also be estimated numericall y: . . . )(2 . . . )(4 . . . (1 (1 )2 log log 2
Page 13
Conservation Physical principles should apply at the discrete level: if ma ss, momentum and energy are conserved, they can only be distributed improp erly Integral form of a generic conservation law ∂t u dV dS q dV, accumulation influx source/sink flux function Caution: nonconservative discretizations may produce reas onably looking results which are totally wrong (e.g. shocks moving with a wro ng speed) even nonconservative schemes can be consistent and stable correct solutions are recovered in the limit of very fine grids Problem: it is usually unclear whether or not the mesh is sufficien tly fine
Page 14
Discrete conservation 1. Any finite volume scheme is conservative by construction both locally (for every single control volume) and globally (for the whole domain) 2. A finite difference scheme is conservative if it can be written in the form +1 +1 which is equivalent to a vertex-centered finite volume discre tization 3. Any finite element scheme is conservative, at least globally =1 ∂u ∂t dV = 0 , i = 1 , . . . , N Summation over yields a discrete counterpart of the integral conservation law ∂u ∂t dV = 0 ∂t dV dS dV
Page 15
Boundedness Convection-dominated / hyperbolic PDEs Pe , Re spurious undershoots and overshoots occur in the vicinity o f steep gradients quantities like densities, temperatures and concentration s become negative the method may become unstable or converge to a wrong weak soluti on 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.2 0.4 0.6 0.8 low-order high-order 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.2 0.4 0.6 0.8 Idea: make sure that important properties of the exact solutio n (monotonicity, positivity, nonincreasing total variation) are inherited by the numerical one
Page 16
Design of nonoscillatory methods Monotone methods ∂u ∂t ∂f ∂x = 0 +1 ) such that ∂H ∂u i, j Then +1 +1 Example. Let ∂u ∂t ∂u ∂x = 0 be discretized by UDS in space, FE in time +1 = 0 , H ) = Derivatives ∂H ∂u = 1 ν, ∂H ∂u ν, where monotone under the CFL condition 1 (cf. stability analysis) Lax-Wendroff theorem: If a monotone consistent and conservative method converges, then it converges to a physically acceptable weak solution
Page 17
Design of nonoscillatory methods Godunov’s theorem: Monotone method are at most first-order accurate Monotonicity-preserving methods (monotone if linear) +1 +1 i, If the initial data is monotone, then so is the solution at all times It is known that the total variation defined as TV ) = ∂u ∂x dx is a nonincreasing function of time for any physically admissible weak solution Total variation diminishing methods (monotone if linear) TV +1 TV where TV ) = Classification monotone TVD monotonicity-preserving
Page 18
Total variation diminishing methods Harten’s theorem: An explicit finite difference scheme of the form +1 ) + +1 +1 is total variation diminishing (TVD) provided that the coe cients satisfy , c +1 , c +1 Semi-discrete problem du dt +1 = 0 conservation form Idea: switch between high- and low-order flux approximations depending on the local smoothness of the solution so as to enforce Harten’s conditions: +1 +1 + +1 +1 +1 where 0 +1 2 is a solution-dependent correction factor (flux limiter)
Page 19
TVD discretization of convective terms Example. Pure convection equation ∂u ∂t ∂u ∂x = 0 , v > , f vu Linear flux approximations +1 vu upwind difference +1 +1 central difference Smoothness indicator +1 Nonlinear TVD flux +1 vu Φ( )( +1 Harten’s coefficients 2 2 + Φ( Φ( , c +1 = 0 Standard flux limiters: Φ( ) = 1+ Van Leer Φ( ) = max min , r }} minmod Φ( ) = max min 1+ MC Φ( ) = max min min , r }} superbee 1D stencil i−1 i+1 i+1 i−1