Bernatz R. Fourier Series and Numerical Methods for Partial Differential Equations

Подождите немного. Документ загружается.

232 THE FINITE DIFFERENCE METHOD

Centered differences will be used for both second partíais in the PDE. The spatial

domain (0<x<l)is divided into N equal subintervals. The length of each

interval is h

—

1/N. The time-independent boundary conditions at x = 0 and x = 1

mean new values for y need to be calculated for y^

— 2, 3, ... N

—

1 and N

only. See Figure 9.5 for details.

¿=123...

N N + 1

—I—I—I—I—I—I—I—I—I—I—I—I

►

x= 0 ... 1

Figure 9.5 Node indexing for the ID wave equation.

An implicit procedure will be used, and the corresponding computational cell for

an arbitrary location is shown in Figure (9.6). The resulting finite difference version

i-1

+ 1

n+l

Vi-1

y'i

y?-

f- n+1

n-1

Figure 9.6 Computational domain for the ID wave equation.

of the wave equation is

M h

)

(9.20)

Each term with superscript "n + 1" is an unknown. The equation above is solved for

these terms, and the result is

,n+l ^„,

o».

OÎ/Γ-Τ + (1 + MvT" - ÛÎ/Γ+Τ =

y? - y

(9.21)

where

a = c

1D WAVE EQUATION 233

Equation (9.21) is constructed for each of the interior spatial nodes, and the resulting

system of equations is represented in matrix form as

1 + 2a -a 0 ·

-a 1 + 2a -a 0

0 -a 1 4- 2a -a 0

where

-a

+ 2a J

92 = 2yZ-y%-

+ay

1 Ί

L i/jv-i

9j = 2y?-j£

j = 3...JV-l

= ZVN - VN

2/iv+i

The coefficient matrix is tridiagonal as in the case of the implicit form of the ID heat

equation.

9.7.2 Initial Conditions

We know that the initial conditions for a ID wave equation include both an initial

displacement y(x, 0) and an initial velocity yt(x, 0). There is more than one way that

these conditions may be translated to a finite difference scheme. One such scheme is

presented in this section.

Referring to Figure 9.7, an arbitrary spatial location with index i is considered.

The time "level" of n = 1 corresponds to t — 0. The value of y\ is determined

directly from the prescribed initial condition of the IB

VP.

That is,

Ι/ί=Λ((<-1)Λ)

2,3,

...AT

The node is "filled" in the figure to indicate it is a known quantity. The values for y\

and y

are determined by the respective boundary conditions.

The value for yf is prescribed next by using the initial velocity given by /^(^)·

The simple forward difference of 0(r) is used to write

which is simplified to give

VÏ = Vi + rf

((i - 1)Λ)

-1)Λ)

:2,3,..JV

234 THE FINITE DIFFERENCE METHOD

i-l

t + 1

Yvt

-Φ-

■

n = 2

-Φ-

Ltf

[iff

■n= 1

n = 0

Figure 9.7 Initial condition for the ID wave equation.

As in the case for n = 1, the values for y\ and y\ are set using the respective

boundary conditions.

Now that both initial conditions have been incorporated in the finite difference

scheme, the semi-implicit formula of Equation (9.21) is used to determine future

values for all interior nodes on time planes corresponding to n > 3.

■ EXAMPLE 9.2

The implicit formula developed in Section 9.7.1 is used here to solve IBVP

(9.22).

(

(x,t)

= 2y

{x,t) (PDE)

0 < x < 10, t > 0

IBVP<

y(*,o) = o (ici)

î/t(:r,0)=0 (IC2)

2/(0,

ί) = ^sint (BC1)

l 2/(10,í) = 0 (BC2)

(9.22)

Figure 9.8 shows the string's displacement for various times t.

9.8 2D HEAT EQUATION IN CARTESIAN COORDINATES

The finite difference formulation for the general 2D heat equation given in IBVP

(9.23) is developed in this section. The formulation will be for the case of a homo-

2D HEAT EQUATION IN CARTESIAN COORDINATES

0.5

y o

-0.5

-1

(a)t=1.0

■(b) i = 5.0

9.0

Figure 9.8 Plots of

the

finite

difference solution y versus t for various t.

geneous PDE, where q(x,y,t)

—

(x, y, t) = fcV

w(x, y, t) + g(x, y, Í)

u(x,2/,0) = /(x,2/)

(PDE)

(IC)

IBVP<

aiw(x, 0, t) + α2%(#,0,ί) = g\{x,t) 0<#<c (BCl)

fritte, y, t) + 6

i¿x(c, y,í) = #

(^¿) 0 < y < d (BC2)

c\u(x,d,t)+ C2Uy{x,d,t) = g${x,t) 0<x<c (BC3)

diw(0,j/,i)+d

wx(0,y,i) = </4(y,i) 0<y<d (BC4)

236 THE FINITE DIFFERENCE METHOD

We begin by considering the explicit finite difference form of the 2D, unsteady

heat equation. Replacing the partíais in the PDE above with centered differences in

space and forward differences in time gives

i ,3 i, 3

M-lj

2?y

i,3

-1,3

\j+l

y-i

(9.24)

where uniform grid spacing in both the x and y directions is assumed. The dependent

variable u is known at all nodal locations in the t = n time plane in Figure (9.9). The

only unknown in Equation (9.24) is u^j

, so this equation may be solved explicitly

for that quantity.

+ l

y y y y y

/ / / / /

/ y y y y y

/ y /

y y y y y' y . \

y y y y y y

J 1

y.. / / / / /

/ i-\ i z+1

i + i

Figure 9.9 Node configuration for the explicit

version of

the

2D heat equation.

The fully implicit form of this difference equation would have all terms in u on

the right-hand side taken from the t

—

n +

time plane. Each would be an unknown,

and the finite difference equation would have a total of five unknown quantities. A

system ofmxn equations would result in the fully implicit case (m is the number

of nodes in the x direction, and n is the number of nodes in the y direction). Not

only would the system be large for large m and n, but the matrix would not be the

simple tridiagonal form as in implicit schemes for the ID case. It would be sparse

(i.e.,

most terms are zero), but still costly to solve.

One compromise scheme is the Alternating Direction Implicit (ADI) method.

Implicit calculations in x— and

y—

directions alternate from one time step to another.

For example, we can calculate new values for u on the t = n 4- 1 time plane using

the equation

n+1

- υ

■

^ffij-K*

«£?-*,,

j+1

-2<¿

%j-i

h? h?

(9.25)

2D HEAT EQUATION IN CARTESIAN COORDINATES 237

which is implicit in x only. This semi-implicit equation has only three unknowns.

The u quantities along a given x line are taken from the t = n time plane, and are

known values. See Figure 9.10 for a visual understanding of the computation cell for

this case. The new values of u are calculated in a single step for the single line using

n-f

/ / / / /

+· x

Figure 9.10 Node configuration for a semi-implicit

version of

the 2D

heat equation.

a system of equations and the tridiagonal matrix as in the ID implicit case. Such

a method has to be used for each line on the new time plane. Therefore, the new

values for the dependent variable on the new time plane are calculated by "sweeping"

through the plane in either the x direction or the y direction. Typically, the sweep

direction alternates from one time step to the next, hence the name "alternating."

This semi-implicit method must use an iterative procedure in order to calculate

new values for u for the entire time plane. Here is why. The first row for which u

is changeable is that directly adjacent to the boundary. We will refer to this row as

"row 2." The boundary row is "row 1" and the second interior row is "row 3," and

so on. Rows 1 and 3 are used to calculate new values for row 2, but must be treated

as known quantities. Well, row 3 is a variable row, so it must be calculated as well.

So,

after row 2 is calculated, we use row 2 to calculate row 3, and row 3 will most

likely be different than it was when we used it to calculate row 2. So the values we

calculated for row 2 are no longer valid because row 3 is now different than it was

when we used it to calculate row 2. So, after we are done sweeping through the

t = n + l time plane, we must go back and sweep again using the new values for u

on adjacent rows to calculate u values on a given row. You can see that this process

must be repeated again and again because it suffers from the same problem on each

sweep. Luckily, this iterative process, sweeping time and time again, converges for

the whole time plane to values of u that solve the semi-implicit equation and the fully

implicit form if it were used.

238 THE FINITE DIFFERENCE METHOD

The iterative process described above cannot continually use u values from the

previous time step (the nth time level) in the iterative process to determine u values in

the n

4-1

time

plane.

Instead, intermediate iterative

values,

denoted by the superscripts

m and

m4-1

in Equation 9.26, are used in the semi-implicit

process.

New intermediate

quantities are denoted with the m +

superscript. Values determined in the previous

iteration within this time step calculation are denoted with the "m" superscript.

-aut\j + (1 + 4a)u™

- au^ij = au&

+ au™^ + u^ (9.26)

The iterative process for calculating u vales on the n 4-1 time plane continues until

\\uZ

-u™\\<e

where e is a small, positive real number serving as a convergence criteria. Once the

convergence criteria has been satisfied in the iterative process, the w^t

= u™^

1 anc

calculation for u values on the next time plane begins with the iterative process and

a line-by-line sweep in the alternate direction.

■ EXAMPLE 9.3

The semi-implicit finite difference formulation derived in Section 9.8 is applied

to a semi-homogenous Dirichlet heat problem in this section. The specifics of

the problem are given in IBVP (9.27).

IBVP^

In order to demonstrate the accuracy of the time-dependent solution, we

begin with the function

u(x,y,t) — 10(x

—

)(y

—

y) sint

which satisfies the given initial condition and boundary conditions in (9.27),

and substitute for u in the PDE. Doing so will determine the source function q.

In this case,

q(x, y, t) = 10(x - x

){y

- y) cost - 10

♦

0.1(2x

+ 2y- 2x

- 2y

) sint.

The finite difference solution for u is plotted in Figure 9.11 using a spherical

symbol for nodal values. The finite difference solution was calculated on a

spatial grid with uniform spacing of

0.05.

The time step size in the calculation

process was

0.1.

The true solution u(x, y, t) is plotted in the wire frame. Both

are surfaces are for time t = 3.1. The finite difference solution surface matches

the true solution quite well, even after

time steps. However, it does lag (i.e.,

is lower than) the true solution surface.

(x,y,t)=0AV<

u(x,y,0) =0

u(x,0,t) = 0

tUl,î/,*) = 0

u(x, l,i) = 0

u(0,y,t) = 0-

u(x

q(x

0<£< 1

0<y < 1

0<x < 1

0< y< 1

(PDE)

(IC)

(BC1)

(BC2)

(BC3)

(BC4)

(9.27)

TWO-DIMENSIONAL WAVE EQUATION

239

-0.005

-0.01

-0.015

-0.02

-0.025

-0.03

0.8

0.6

0.4

0.2

0.4

0.8

Figure 9.11 Finite difference solution (point plot) and true solution (wire frame).

9.9 TWO-DIMENSIONAL WAVE EQUATION

The general IBVP

for a

wave on

2D elastic surface stretched from

0 < x < L

and

< y < M

is given in IBVP (9.28).

this section,

semi-implicit finite difference

formula

for

the nonhomogeneous PDE

the IBVP

derived.

IBVP<

(x,y,t)--

< x < L,

--c

z(x,y,t)+F(x,y,t)

< y < M, t > 0

z(x,y,0)

fi(x,y)

(x,y,0)

(x,y)

aiZy(x,0,t)

biz(x,0,t)

gi(x,t)

(L,

y, t) +

z(L,

y, t) = g

(y, t)

Zy(x, M,

t) +

z(x,

t) = g

(x, t)

^4^(0,

y, t) +

z(0,

y, t) = g

(y, t)

(PDE)

(ICI)

(IC2)

(BC1)

(BC2)

(BC3)

(BC4)

(9.28)

Figure

9.12

shows

the

three time planes involved

approximating

the

second

partial

of z

with respect

time

t in

the semi-implicit formulation. The values

of z

that are to be determined and are considered as unknown are those on the

+1 plane.

The resulting finite difference version

the homogeneous PDE is given in Equation

(9.29).

+ z

n+l

„n+1

-2z"

n+l

i+lj

-1

ij ^

lj+l

(9.29)

The three unknown quantities, those with superscripts

of n

left-hand side. All other terms are moved

the right-hand side. Letting

(assuming

h = Ay),

Equation (9.29) simplifies

are moved to the

az^. +

2a)zZ

n+l

i+ij

α4_

(2-2α)4 +

α4

-4-

(9.30)

9.10 2D HEAT EQUATION IN POLAR COORDINATES

The finite difference solution

for

the heat equation

polar coordinates

presented

in this section. The resulting formulation uses the semi-implicit ADI method where

240 THE FINITE DIFFERENCE METHOD

n+1

n-1

//////

2—1

2 2+1

Figure 9.12 Cell configuration for

semi-implicit

version of

the

2D wave equation.

unknown dependent variable values

are

calculated either along

line

constant

angle

or constant radius

The equation

for

heat u(r, 0, t) in polar coordinates

given

Equation (9.31).

2 βφ2

~dt

u 1

r δφ r*

q(r,<l>,t)

(9.31)

The thermal diffusivity is represented by

and the source term is given by q(r, φ,

t).

Figure 9.13 depicts an arbitrary finite difference element in polar coordinates with

uniform grid spacing

both

and

φ.

The constant radius contours are indexed

and the

constant

contours

by i The

superscript

denotes

the

corresponding

value of u is taken from the

t = n

time plane, while the superscript

n +

implies the

corresponding value

of u

is an unknown value

for u at

that location in the

t = n + 1

time plane.

The

formulation

for the

finite difference expression corresponding

Equation (9.31) will be for semi-implicit 'updating" for

lines

of constant φ

indicated

in the FD element

Figure 9.13.

Replacing

the

partial derivatives

Equation (9.31) with finite difference repre-

sentations gives

,.η+ι

i>>3

n+l

n+l i

n+l

.n+1

%j+l

■Ui

n+1

j-l

(Ar)

jAr

(Αφ)

2Ar

€

i,3

(9.32)

Note: Each spatial finite difference expression in Equation (9.32) has truncation error

of second order, either

in Ar or Αφ.

The temporal derivative has

error

order

2D HEAT EQUATION IN POLAR COORDINATES 241

.i + i

,n+l

'»+1.J

— i ¿ ¿4-1

Figure 9.13 Finite difference element for polar coordinates.

Multiplying both sides of Equation (9.32) by At and rearranging terms u by

coordinate indices and time plane, the result is

-CSU^+CU^-CNU^ = CwuUj + Cmuli

+hj

+ Atq^j (9.33)

where expressions for the coefficients Cs, CN, ... are given in Table 9.1.

Table 9.1 Finite Difference Coefficients for Equation (9.33)

where

Atk

(Ar)

C = 2a + 1

c -i

2ß

and

—

OÍ

4- —

¿3

β-

* (Δ0)

A similar semi-implicit finite difference representation for Equation (9.31) can

be derived for the case where updating is done along contours of constant r. The

simplified algebraic equation is

-CWU^IJ + CU^-CEU^ = Csu^+Cmuîj

The only coefficient formulas for Equation (9.34) different from those in Equation

(9.33) are those for C and C

, which become

C = 1 + -hr and C

= 1 - 2a

There are two aspects of the polar coordinate finite difference formulation that

need attention. The first concerns the finite difference solution of Equation (9.31) at