Loading...

PPT – Numerical Methods for Partial Differential Equations PowerPoint presentation | free to download - id: 550d95-MWNlY

The Adobe Flash plugin is needed to view this content

Numerical Methods for Partial Differential

Equations

- CAAM 452
- Spring 2005
- Lecture 7
- Instructor Tim Warburton

Summary of Small Theta Analysis

- The dominant remainder term in this analysis

relates to a commonly used, physically motivated

description of the shortfall of the method

Dissipative Unstable Dispersive Dispersi

ve

Dissipation in Action

- Consider the right difference version
- We are going to experimentally determine how much

dissipation the solution experiences.

Testing Methodology

- Reduce the time-stepping error to a secondary

effect by choosing a 4th order Runge-Kutta (JST)

time-integrator and a small dt. - Fix the method (choose one of the difference

operators for the spatial derivative) - Do a parameter study in M, i.e. we ask the

questions how does increasing the number of data

points change the error. - We need to understand what questions we are

asking - Is the computer code stable (as predicted by the

theory)? - Does the computer code converge (as predicted by

the theory)? - What order of accuracy in M do we achieve?
- We hypothesize that if the theory holds then we

should achieve - What is the actual approximate p achieved?

Initial Condition

- Because we do not wish to introduce uncertainty

over the source of errors in the computation we

use an initial condition which is infinitely

smooth.

Computing Approximate Order of Accuracy

- We compute the error by
- We will compute the approximate order of accuracy

by assuming - where C is independent of M
- If we measure the error for two different M then

we can eliminate the constant - In the case when M22M1

After 4 Periods

- The numerical pulses are in the right place but

have severely diminished amplitude.

After 4 Periods

M Maximum error at t8 order

20 0.69149169495179

40 0.69137534363447 0

80 0.68329408124856 0.01

160 0.63163960492185 0.11

320 0.52742560650408 0.26

After 4 periods the solution is totally flattened

in all but the last 2 results. If we botheredto

keep increasing M we would eventually see the

error decline as 1/M

The Unstable Left Stencil

- I repeated this with everything the same, except

the choice of instead of

We clearly see that there is initial growth near

the pulses, but eventually the dominant feature

is the highly oscillatory and explosive growth

(large m in the above red term).

Cont (snapshot)

Dispersion in Action

- Consider the central difference version
- With the same time integrator before, M100,
- We note that the remainder terms are dispersive

corrections i.e. they indicate that modes of

different frequency will travel with different

speed. - Furthermore, to leading order accuracy the higher

order modes will travel more and more slowly as m

increases.

2nd Order Central DifferenceAfter 4 Periods

- We notice that the humps start to shed rear

oscillations as the higher frequency Fourier

components lag behind the lower frequency Fourier

components.

Convergence Study (t8)

What should we use as an error Indicator ??

2nd Order Central Difference

M Maximum error at t8 order

20 0.94240801474506

40 0.35029899850733 1.42

80 0.38298512776859 0.13

160 0.19118127466878 1.00

320 0.05471730524482 1.80

- We do not see a convincing 2nd order accuracy
- I computed this by log2(error_M/errorM-1)

4th Order Central Differencing

4th Order Central Difference

M Maximum error at t8 order

20 0.44353878542277

40 0.12852952972116 1.79

80 0.02040190886767 2.66

160 0.00134076478259 3.93

320 0.00008319438414 4.01

- We see pretty convincing 4th order accuracy
- I computed this by log2(error_M/errorM-1)

6th Order Central Difference

M Maximum error at t8 order

20 0.19409575615520

40 0.04953652526983 1.97

80 0.00169299155603 4.87

160 0.00002891724105 5.87

320 0.00000046041415 5.97

- We see pretty convincing 6th order accuracy
- I computed this by log2(error_M/errorM-1)

Summary of Testing Procedure

- Understand what you want to test
- Keep as many parameters fixed as possible
- If possible, perform an analysis before hand
- Run parameter sweeps to determine if code agrees

with analysis. - Estimate error scaling with a single parameter if

possible. - It should be apparent here that the errors

computed in the simulations only approximate the

tidy power law (dxp) in the limit of small dx

(large M). - It is quite common to refer to the range of M

before the error scaling applied as the

preasymptotic convergence range. - The preasymptotic range is due to the other

factors in the error estimate which are much

larger than the dxp term (i.e. the pth

derivatives may be very large, or the constant

may be large) - In this case we heuristically set dt small, using

a high-order Runge-Kutta time integrator, and

then performed a parameter sweep on M. We could

have been unlucky and the time-stepping errors

may have been dominant which would mask the dx

behavior. - To be careful we would try this experiment with

even smaller dt and check if the scalings still

hold.

Summary of our Approach to Designing Finite

Difference Methods

- We have systematically created finite difference

methods by separating the treatment of space and

time derivatives. - Then designing a solver consists of

choosing/testing - a time integrator (so far we covered

Euler-Forward, LeapFrog, Adams-Bashford(2,3,4),

Runge-Kutta) - a discretization for spatial derivatives
- a discrete differential operator which has all

eigenvalues in the left half of the complex plane

(assuming the PDE only admits non-growing

solutions). - Possibly using Gerschgorins theorem to localize

the eigenvalues of the discrete differentiation

operator - dt so that dtlargest eigenvalue (by magnitude)

of the derivative operator sits inside the

stability region ? stability. - (small dx) spatial truncation analysis ?

consistency and accuracy. - Fourier analysis for classification of the

differential operator. - Writing code and testing.

Some Classical Finite Difference Methods

- However, there are a number of classical methods

which we have not discussed and do not quite fit

into this framework. - We include these for completeness..

Leap Frog (space and time)

- We use centered differencing for both space and

time. - We know that the leap frog time stepping method

is only stable for operators with eigenvalues in

the range - However, we also know that the centered

difference derivative matrix is a skew symmetric

matrix with eigenvalues - So we are left with a condition

i.e.

cont

- We can perform a full truncation analysis (in

space and time) - We know that dt lt dx so

Lax-Friedrichs Method

- We immediately conclude that the following method

is not stable - Because the Euler-Forward time integrator only

admits one stable point (the origin) on the

imaginary axis, but the central differencing

matrix has all eigenvalues on the imaginary axis.

- However, we can stabilize this formula by

replacing the second term in the time-stepping

formula

cont

- This formula does not quite fit into our

constructive process (method of lines approach). - We have admitted spatial averaging into the

discretization of the time derivative. - We can rearrange this

cont

- We can immediately determine that this is a

stable method as long as cdt/dx lt1 - Given, this condition we observe that the time

updated solution is always bounded between the

values of the left and right neighbor at the

previous time because this is an interpolation

formula.

cont

- A formal stability analysis would involve
- Which gives stability for each mode if

cont

- Thus considering all the possible modes
- We note that the middle mode requires
- This condition gives a sufficient condition for

all modes to be bounded. - By the invertibility and boundedness of the

Fourier transform we conclude that the original

equation is stable.

cont

- We can recast the Lax-Friedrichs method again
- This method consists of Euler Forward, central

differencing for the space derivative and an

extra dissipative term (i.e. a grid dependent

advection diffusion equation

)

CFL Number

- The ratio appears repeatedly, in

particular in the estimates for the maximum

possible time step. - We refer to this quantity as the CFL

(Courant-Friedrichs-Lewy) number. - Bounds of the form which result

from a stability analysis are frequently referred

to as CFL conditions.

Lax-Wendroff Method

- We are fairly free to choose the parameter in the

stabilizing term - The artificial viscosity term acts to shift the

eigenvalues of the spatial operator into the left

half-plane. - Recall the Euler-Forward stability region is

the unit circle centered at -1 in the complex

plane. So pushing the eigenvalues off the

imaginary axis allows us to choose a dt small

enough to push the eigenvalues of the discrete

spatial operator into the stability region

Class Exercise

- Please provide a CFL condition for
- In terms of
- What is the truncation error?

Von Neumann Analysis

- By stealth I have introduced the classical Von

Neumann analysis. - The first step is to Fourier transform the finite

difference equation in space. - A short cut is to make the following

substitutions

Analysis cont

- We can also make the substitutions for the

difference operators

Example

- The Lax-Wendroff scheme
- becomes

Example cont

- So the Lax-Wendroff scheme becomes a neat one

step method for each Fourier coefficient - We have neatly decoupled the Fourier coefficients

so we are back to solving a recurrence relation

in n for each m - Hence we need to examine the roots of the

stability polynomial for the above

Example cont

- This case is trivial as we just need to bound the

single root by 1 - Plot it this as a function of theta

Final Von Neuman Shortcut

- We can skip lots of steps by making the direct

substitution - Here g is referred to as the amplification factor

and imtheta is our previous theta_m

Next Time

- Definition of consistency
- Lax-Richtmyer equivalence theory
- Treating higher order derivatives

Homework 3

- Q1) Build a finite-difference solver for
- Q1a) use your Cash-Karp Runge-Kutta time

integrator from HW2 for time stepping - Q1b) use the 4th order central difference in

space (periodic domain) - Q1c) perform a stability analysis for the

time-stepping (based on visual inspection of the

CK R-K stability region containing the imaginary

axis) - Q1d) bound the spectral radius of the spatial

operator - Q1e) choose a dt well in the stability region
- Q1f) perform four runs with initial condition
- (use M20,40,80,160) and compute maximum

error at t8 - Q1g) estimate the accuracy order of the solution.
- Q1h) extra credit perform adaptive time-stepping

to keep the local truncation error from time

stepping bounded by a tolerance.