Unformatted text preview: .3.3. The LaxWendro scheme for the kinematic wave equation ut + aux = 0
is obtained from the Taylor's series in time
2
un+1 = un + t(ut)n + 2t (utt )n + : : : :
j
j
j
j Time derivatives are replaced by spatial derivatives that are obtained from the partial
di erential equation. For simplicity, let us suppose that a is a constant, then ut = ;aux utt = ;auxt = ;a(ut )x = a2uxx ::: 3.3. Matrix Stability Analysis 25 (In fact, the LaxWendro scheme is rather di cult to implement when is not constant.
We'll discuss an alternative in Chapter 6 that is superior to this approach in this case.)
The Taylor's series is truncated and the spatial derivatives are approximated by centered di erences (2.1.7, 2.1.9). Retaining O( t2 ) terms leads to the secondorder LaxWendro scheme
2 or Ujn+1 = Ujn ; 2 (Ujn+1 ; Ujn;1) + 2 (Ujn+1 ; 2Ujn + Ujn;1)
Ujn+1 = ( 2+ 1) Ujn;1 + (1 ; 2)Ujn + ( 2; 1) Ujn+1: (3.3.8) where, as usual, is the Courant number. Taylor's series arguments immediately show
that the local error is O( t2) + O( x2 ).
The coe cients of this LaxWendro scheme are not all positive. For example, if
0 < < 1 then the coe cient of Ujn+1 will be negative. Likewise, the coe cient of Ujn;1
is negative when is negative. Thus, Theorem 3.1.1 cannot be used to establish stability
in the maximum norm. As expected, the LaxWendro scheme gives the exact solution
of the kinematic wave equation when j j = 1.
We'll establish the stability of a periodic initial value problem in the Euclidean norm
when j j 1 using de nition (3.1.14b). The von Neumann method could also be used
in this case (cf. Problem 3.2.3). To begin, we square (3.3.8) to obtain
(Ujn+1 )2 = ( + 1) Ujn;1 + (1 ; 2 )Ujn + ( ; 1) Ujn+1]2:
2
2
If j j 1 then
2 (1 ; 2 )
Ujn;1 ; 2Ujn + Ujn+1]2 0:
4
Add this term to the right side of the previous expression and sum over a period to get
J ;1
J ;1
2
X n+1 2 X 2(1 + ) n 2
(Uj;1) + (1 ; 2)(Ujn )2 + (12; ) (Ujn+1)2
(Uj )
2
j =0
j =0 ; (1 ; 2)(Ujn+1Ujn ; Ujn Ujn;1)]: This expression can be simpli ed by reindexing the summations, e.g.,
J ;1
X
j =0 (U ;1
n
j )2 = J ;2
X k (U =;1 n
k J ;1
2 = X(U n )2 + (U n )2 ; (U n )2 :
)
k
;1
J ;1
k =0 26 Basic Theoretical Concepts n
n
The solution is periodic so U;1 = UJ ;1 hence,
J ;1
X
j =0 (Ujn;1)2 = Similar reindexing of other terms yields
J ;1
X J ;1
2 = X(U n )2
(U +1)
k
j =0
k =0
n
j J ;1
X
k =0 J ;1
X
j =0 (Ukn)2 : U U ;1 =
n
j n
j J ;1
X
k =0 Ukn+1Ukn : Thus, the summation simpli es to
J ;1
X
j =0 (U n
j +1 )2 J ;1
J ;1
2
X 2(1 + )
2 ) + (1 ; ) ](U n )2 = X(U n )2 :
+ (1 ;
j
j
2
2
j =0 j =0 Therefore, kUn+1 k2 kUnk2 and the LaxWendro scheme is stable when the Courant,
Friedrichs, Lewy Theorem (j j 1) is satis ed
Centered schemes like the LaxWendro (3.3.8) and LaxFriedrichs (3.2.7) methods
may require arti cial boundary conditions for initialboundary value problems. For instance, suppose the kinematic wave equation with a positive wave speed a is to be solved
on 0 < x < 1. This problem only requires a boundary condition at x = 0 (Section 1.3). A
numerical solution at the rightmost point j = J , however, cannot be computed by either
n
the LaxWendro or LaxFriedrichs schemes. Another method is needed to compute UJ .
As we'll learn in Chapter 6, a poor choice could a ect the stability of the method. Problems 1. Write a computer program for the di erence scheme (3.2.8). Implement it with
= 2, which corresponds to the LaxWendro scheme. Assume that the initial
data
u(x 0) = (x)
is periodic in x with period 2. The problem should be solved on ;1
0 < t T . Use J (= 2= x), , and N (= T= t) as input parameters.
1.1. Execute your program when a = 1 and (x) has the form
8
< x + x if ; x x < 0
(x) = : x ; x if 0 x < x
:
0
elsewhere for x 2 (;1 1) x 1, 3.4. The Lax Equivalence Theorem 27 1.2. This data is an attempt to simulate the e ects of a small error introduced by,
say, round o . Execute your program for about 20 time steps with J = 10 and
= 0.5, 0.999, 1.1 (more if you like). Plot the numerical and exact solutions
as functions of x for a few times. Comment on the solutions for each value
of . Which choice of would be preferred if the initial data actually did
correspond to a rounding error?
1.3. Solve a problem...
View
Full Document
 Spring '14
 JosephE.Flaherty
 Numerical Analysis, Fourier Series, Partial differential equation, Hilbert space, di erence, nite di erence

Click to edit the document details