Download Explicit Euler Method for Solving the Diffusion Equation: Stability Analysis and more Study notes Calculus in PDF only on Docsity!
Diffusion Equation and Final Exam Review
Adrian Down
May 04, 2006
1 Implicit Euler solutions to the diffusion equa-
tion
1.1 Review
Last time, we developed a numerical approach to solve the diffusion equation,
ct − Dcxx = 0
The method approximated the derivatives using finite difference approxima- tions in analogy to the Forward Euler method. The points used to calculate the difference approximations were on a mesh that had regular spacing ∆x in the space component and ∆t in the time component. We saw that this method become divergent if the spatial mesh spacing was made too small. In particular, we found a condition for stability,
D∆t (∆x)^2
The stability of this method is not ideal. To develop a more stable method, we consider a different finite difference approximation to the deriva- tives.
1.2 Explicit Euler approximation method
1.2.1 Motivation
Our previous method of approximating the finite difference approximations used three spatial mesh points at a particular time, (x − ∆x, t), (x, t), and
(x + ∆x, t), to compute the value of the approximation at one spatial point at the next time step, (x, t + ∆t). Each set of three spatial mesh points contributes to the next time step at one spatial point. Instead, consider taking only one mesh point, (x, t), and from it calcu- lating the point at the next time step, x, t + ∆t), as well as the neighboring points (x − ∆x, t + ∆t) and (x + ∆x, t + ∆t). (x, t∆t) is related to the neigh- boring spatial mesh points by the difference approximation, which defines a system that must be solved for (x, t+∆t) knowing (x, t). In matrix form, the resulting system is tri-diagonal, which can be solved using methods already discussed in this course.
1.2.2 Calculation
Suppose we have the points γm,n, from which we would like to calculate γm− 1 ,n+1, γm,n+1, and γm+1,n+1. The spatial mesh points at t + ∆t are related by the difference approximation,
cxx ≈
γm+1,n+1 − 2 γm,n+1 + γm− 1 ,n+ (∆x)^2
The points γm,n and γm,n+1 are related in the usual way by the finite difference approximation to the time derivative,
ct ≈
γm,n+1 − γm,n ∆t
Setting these equations equal defines the approximation,
γm,n+1 − γm,n ∆t
D
(∆x)^2
{γm+1,n − 2 γm,n + γm− 1 ,n}
⇒ γm,n+1 − γm,n =
D∆t (∆x)^2
{γm+1,n − 2 γm,n + γm− 1 ,n}
The crucial parameter for stability analysis is the prefactor,
D∆t (∆x)^2
Taking k = (^) ∆πx , the sine term in the dispersion relation equals 1, and the dispersion relation becomes,
1 − e−σ∆t^ = −
4 D∆t (∆x)^2
Isolating σ∆t,
e−σ∆t^ = 1 +
4 D∆t (∆x)^2
≡ μ
As before, the quantity on the right, which we define as μ, determines the stability of the method.
Note. In this case, μ is positive definite, in particular, μ > 1, in contrast to the case of the Explicit Euler method, where μ could take negative values. This sign difference will be crucial for ensuring the stability of the Explicit Euler method.
1.3.3 Complex σ∆t
As before, we examine the stability for complex values of σ∆t. If we can show that the real part of σ∆t is less than 0, then any disturbances in the plane wave solution will decay exponentially, giving a stable method. Otherwise, perturbations will grow without limit, and the method will not be stable. Write σ∆t in terms of real and imaginary components, σ∆t = S + ıW. Substituting into the dispersion relation,
e−(S+ıW^ )^ = e−S^ cos W − ıe−S^ sin W = μ
Match the real and imaginary components of the above equation, noting that μ is purely real,
e−S^ cos W = μ e−S^ sin W = 0
The second equation determines W. Since the exponential factor cannot be 0 for any value of S, sin W must equal 0,
sin W = 0 ⇒ W = nπ ⇒ cos W = (−1)n ⇒ e−S^ (−1)n^ = μ > 1
If n is odd, there is no solution for S, and hence such modes are not possible. In the case of n even,
e−S^ = μ > 1
Since ex^ > 1 ⇔ x > 0, it must be that −S > 0 ⇒ S < 0. Hence the plane wave solutions decay exponentially, and the method is stable.
Note. We have just shown that the Implicit Euler method is stable for all values of (^) (∆D∆x)t 2 , whereas we saw that the Explicit Euler method was only
stable for (^) (∆D∆x)t 2 < 12. Hence the Explicit Euler method is in general more stable, although it involves more computation to solve the tri-diagonal system for each iteration.
2 Final exam review
2.1 Space-time coordinates
Time 8:00 AM, May 18
Location 4 LeConte
2.2 General overview of final topics
Problem 1: Richardson extrapolation This method is used to create successively more accurate approximations by taking linear combina- tions of approximations with different mesh sizes. Be familiar with the following:
- Recognize when this method can be used. Richardson extrapo- lation can be used to approximate any function that can be ex- pressed as a Taylor series in only even powers of a variable.
- How to construct a Richardson table.
Problem 2: Trapezoidal and Simpson’s rules Be familiar with the fol- lowing:
- Simpson’s rule is an extension of Trapezoidal method.
- How to compute the truncation error of both methods.
Problem 6: PDE difference scheme This problem is given in advance. You can prepare solutions using whatever resources you like. The prob- lem will be to analyze the stability of wave solutions to a particular PDE, the advection equation,
ct + cx = 0
2.3 Example for problem 1: Richardson extrapolation
We see how Richardson extrapolation method can be used to obtain a very accurate approximation to a function. Recall the following exercise from the beginning of the course, in which we compute ln 2 using the derivative of 2x^ evaluated at x = 0,
ln 2 = (2x)′(0)
We could write the derivative as a centered difference quotient,
ln 2 = lim h→ 0
2 h^ − 2 −h 2 h
We can create a function by defining a value at the removable singularity at h = 0,
f (h) =
2 h− 2 −h 2 h h^6 = 0 ln 2 h = 0
Notice that f (h) is an even function, and hence its Taylor series contains only even powers of h. This information is sufficient to tell us that Richardson extrapolation can be used to approximate the function; we do not need to know the exact coefficients of the expansion. Build the Richardson table using the algorithm developed previously. Each successive level of the Richardson table cancels a further term in the error of the approximation,
f 1 (h) =
f
h 2
f (h) = ln 2 + O(h^4 )
This linear combination has removed the second order error. We can continue further to remove higher order terms,
f 2 (h) =
f 1
h 2
f 1 (h) = ln 2 + O(h^6 )
Computing f 1
(h 2
implicitly requires us to evaluate f
(h 4
. A table format is useful to evaluating all of the necessary function values. This table is known as a Richardson table. As an example, consider a step size of h = 1, which is a very poor approx- imation to taking the limit as h → 0. We will see, however, that even with this poor step size, the approximation can be made very accurate. Taking
f (1) = 34 =. 750000
− (^13) −−→ f 1 (1) =. 692809
− 151 −−→ f 2 (1) =. 693147 (^43) ↗
(^1615) ↗ f
2
= √^12 =. 707107
− (^13) −−→ f 1
2
4 3 ↗ f
4
Table 1: Richardson table for f (h) ≈ (2x)′(0)
the difference with the exact result,
ln 2 =. 693147 ⇒ f 2 (1) − ln 2 = 2. 4042 × 10 −^7
Hence the result of the Richardson extrapolation is accurate to one part in a million in only three iterations, even though we used a very poor initial mesh size.
2.4 Problem 6: Advection equation
2.4.1 Background
Consider the advection equation,
ct + cx = 0
This equation describes the transport of a material by some type of flow, as opposed to the diffusion equation, which described transport by diffusion. This equation is useful because it allows us to test complex numerical approximation methods with a simple benchmark example. For a numerical approximation scheme to be useful, it should at least be able to solve this simple example.
2.4.4 Stability analysis
You will be given one of these methods, for which you must analyze the stability of the plane wave solutions. The general method is:
- Assume wave solutions of the form,
γm,n = eσ(n∆t)+ık(m∆x)
- Substitute the assumed form of the solution into the approximation scheme.
- Derive the discrete dispersion relation.
- Consider the largest possible perturbations by taking k = (^) ∆^2 πx.
- Write σ∆t in terms of real and imaginary components, σ∆t = S + ıW. Match real an imaginary components of the resulting dispersion relation to solve for W and S.
- Stability is determined be whether S = <(σ∆t) > 0. If S > 0, small perturbations in the initial conditions grow exponentially, and the ap- proximation method is not stable. If S < 0, oscillations decay expo- nentially, and the method is stable.
