416
HYPERBOLIC PDE
PARTIAL DIFFERENTIAL EQUATIONS
417
0.3 0.2 0.2 0 0.2 1 0.5 1 0 0 (a) 0 0.2 0.3 2 0.2 0.4 0.6 (b) A snap shot 0.8 x 1 0.1 0 0.1
We need the solution of this equation to be inside the unit circle for stability, which requires 1 t2 r , r =A 2 1 (9.3.7) 1 cos(/P ) x
The objective of the following MATLAB routine wave() is to implement this algorithm for solving a one-dimensional wave equation.
t
Example 9.4. A Hyperbolic PDE: One-Dimensional Wave (Vibration). Consider a one-dimensional hyperbolic PDE (E9.4.1)
9.3.2 Two-Dimensional Hyperbolic PDE
u(x, t) u (x, t) = x 2 t 2
for 0 x 2, 0 y 2, and 0 t 2
Figure 9.5 A solution for a 1-D hyperbolic PDE obtained by using wave() (Example 9.4).
with the initial conditions and boundary conditions for t = 0 for x = 1 A for (E9.4.2b) (E9.4.2a)
u(x, 0) = x(1 x),
u/t (x, 0) = 0
In this section, we consider a two-dimensional wave equation for the amplitude function u(x, y, t) ((x, y ) is position, t is time) as 2 u(x, y, t) 2 u(x, y, t) + x 2 y 2 = 2 u(x, t) t 2 (9.3.8)
u(0, t) = 0
for x = 0,
u(1, t) = 0
We made the following MATLAB program solve_wave.m in order to use the routine wave() to solve this equation and ran this program to get the result shown in Fig. 9.5 and see a dynamic picture.
0 x xf , 0 y yf , 0 t T
In order for this equation to be solvable, we should be provided with the boundary conditions u(0, y, t) = bx 0 (y, t), u(x, 0, t) = by 0 (x, t), u(xf , y, t) = bxf (y, t), and u(x, yf , t) = byf (x, t)
as well as the initial condition u(x, y, 0) = i0 (x, y) and u/t |t =0 (x, y, 0) = i0 (x, y). In the same way as with the one-dimensional case, we replace the second derivatives on both sides by their three-point central difference approximation (5.3.1) as A
k k uk i,j +1 2ui,j + ui,j 1
x2 with
k k uk i +1,j 2ui,j + ui 1,j
y2 xf x= , Mx yf y= , Ny which leads to the explicit central difference method:
k +1 k 1 ui,j 2uk i,j + ui,j
t2 (9.3.9) T t= N
function [u,x,t] = wave(a,xf,T,it0,i1t0,bx0,bxf,M,N) %solve a u_xx = u_tt for 0<=x<=xf, 0<=t<=T % Initial Condition: u(x,0) = it0(x), u_t(x,0) = i1t0(x) % Boundary Condition: u(0,t)= bx0(t), u(xf,t) = bxf(t) % M = # of subintervals along x axis % N = # of subintervals along t axis dx = xf/M; x = [0:M]*dx; dt = T/N; t = [0:N]*dt; for i = 1:M + 1, u(i,1) = it0(x(i)); end for k = 1:N + 1 u([1 M + 1],k) = [bx0(t(k)); bxf(t(k))]; end r = a*(dt/dx)^ 2; r1 = r/2; r2 = 2*(1 - r); u(2:M,2) = r1*u(1:M - 1,1) + (1 - r)*u(2:M,1) + r1*u(3:M + 1,1) ... + dt*i1t0(x(2:M)); %Eq.(9.3.6) for k = 3:N + 1 u(2:M,k) = r*u(1:M - 1,k - 1) + r2*u(2:M,k-1) + r*u(3:M + 1,k - 1)... - u(2:M,k - 2); %Eq.(9.3.3) end
%solve_wave a = 1; it0 = inline(x.*(1-x),x); i1t0 = inline(0); %(E9.4.2a) bx0t = inline(0); bxft = inline(0); %(E9.4.2b) xf = 1; M = 20; T = 2; N = 50; [u,x,t] = wave(a,xf,T,it0,i1t0,bx0t,bxft,M,N); figure(1), clf mesh(t,x,u) figure(2), clf for n = 1:N %dynamic picture plot(x,u(:,n)), axis([0 xf -0.3 0.3]), pause(0.2) end
k +1 k 1 k k k k = rx (uk ui,j i,j +1 + ui,j 1 ) + 2(1 rx ry )ui,j + ry (ui +1,j + ui 1,j ) ui,j (9.3.10) t2 t2 with rx = A 2 , ry = A 2 x y