Page 429 - Applied Numerical Methods Using MATLAB

418 PARTIAL DIFFERENTIAL EQUATIONS
Since u −1 = u(x j ,y i , − t) is not given, we cannot get u 1 directly from this
i,j i,j
formula (9.3.10) with k = 0:

1 0 0 0 0 0 −1
u i,j = r x (u i,j+1 + u i,j−1 ) + 2(1 − r x − r y )u i,j + r y (u i+1,j + u i−1,j ) − u i,j
(9.3.11)
Therefore, we approximate the initial condition on the derivative by the central
difference as
u 1 − u −1
i,j i,j
= i (x j ,y i ) (9.3.12)
0
2 t
−1
and make use of this to remove u from Eq. (9.3.11) to have
i,j
1
u 1 = {r x (u 0 + u 0 ) + r y (u 0 + u 0 )}
i,j 2 i,j+1 i,j−1 i+1,j i−1,j

+ 2(1 − r x − r y )u 0 + i (x j ,y i ) t (9.3.13)
i,j 0
We use this Eq. (9.3.13) together with the initial conditions to get u 1 and then
i,j
go on using Eq. (9.3.10) for k = 1, 2,. ... A sufficient condition for stability [S-1,
Section 9.6] is
4A t 2
r = ≤ 1 (9.3.14)
2
x + y 2
The objective of the MATLAB routine “wave2()” is to implement this algo-
rithm for solving a two-dimensional wave equation.

Example 9.5. A Hyperbolic PDE: Two-Dimensional Wave (Vibration) Over a
Square Membrane. Consider a two-dimensional hyperbolic PDE

2
2
2
1 ∂ u(x, y, t) ∂ u(x, y, t) ∂u (x, y, t)
+ =
4 ∂x 2 ∂y 2 ∂t 2
for 0 ≤ x ≤ 2, 0 ≤ y ≤ 2 and 0 ≤ t ≤ 2 (E9.5.1)
with the zero boundary conditions and the initial conditions

u(0,y,t) = 0, u(2,y,t) = 0, u(x, 0,t) = 0, u(x, 2,t) = 0 (E9.5.2)
u(x, y, 0) = 0.1sin(πx) sin(πy/2), ∂u/∂t(x, y, 0) = 0for t = 0 (E9.5.3)

We made the following MATLAB program “solve_wave2.m”in order to use
the routine “wave2()” for solving this equation and ran this program to get the
result shown in Fig. 9.6 and see a dynamic picture. Note that we can be sure of
stability, since we have

4A t 2 4(1/4)(2/20) 2 1
r = = = ≤ 1 (E9.5.4)
2
2
x + y 2 (2/20) + (2/20) 2 2