Numerical Methods for Partial Differential Equations - PowerPoint PPT Presentation


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

Get the plugin now

View by Category
About This Presentation

Numerical Methods for Partial Differential Equations


Numerical Methods for Partial Differential Equations CAAM 452 Spring 2005 Lecture 7 Instructor: Tim Warburton Summary of Small Theta Analysis The dominant remainder ... – PowerPoint PPT presentation

Number of Views:204
Avg rating:3.0/5.0
Slides: 40
Provided by: TimWar
Learn more at:


Write a Comment
User Comments (0)
Transcript and Presenter's Notes

Title: Numerical Methods for Partial Differential Equations

Numerical Methods for Partial Differential
  • 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
Dissipation in Action
  • Consider the right difference version
  • We are going to experimentally determine how much
    dissipation the solution experiences.

Testing Methodology
  1. Reduce the time-stepping error to a secondary
    effect by choosing a 4th order Runge-Kutta (JST)
    time-integrator and a small dt.
  2. Fix the method (choose one of the difference
    operators for the spatial derivative)
  3. Do a parameter study in M, i.e. we ask the
    questions how does increasing the number of data
    points change the error.
  4. We need to understand what questions we are
  5. Is the computer code stable (as predicted by the
  6. Does the computer code converge (as predicted by
    the theory)?
  7. What order of accuracy in M do we achieve?
  8. We hypothesize that if the theory holds then we
    should achieve
  9. 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

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
  • Furthermore, to leading order accuracy the higher
    order modes will travel more and more slowly as m

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

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
  1. Understand what you want to test
  2. Keep as many parameters fixed as possible
  3. If possible, perform an analysis before hand
  4. Run parameter sweeps to determine if code agrees
    with analysis.
  5. Estimate error scaling with a single parameter if
  6. 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).
  7. It is quite common to refer to the range of M
    before the error scaling applied as the
    preasymptotic convergence range.
  8. 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)
  9. 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
  10. To be careful we would try this experiment with
    even smaller dt and check if the scalings still

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
  • a time integrator (so far we covered
    Euler-Forward, LeapFrog, Adams-Bashford(2,3,4),
  • 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
  • Possibly using Gerschgorins theorem to localize
    the eigenvalues of the discrete differentiation
  • 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
  • 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

  • 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

  • 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

  • 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

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

  • 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.

  • 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
  • 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

Analysis cont
  • We can also make the substitutions for the
    difference operators

  • 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
  • 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
  • Q1d) bound the spectral radius of the spatial
  • 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.