Partial Differential Equations Notes by Robert Pich´e, Tampere University of Technology

10

Solving Two-Dimensional Laplace Equations • Laplace equation boundary value problems in a disk, a rectangle, a wedge, and in a region outside a circle

10.1

Dirichlet Problem in a disk

Consider the two dimensional Laplace equation in the disk x2 + y 2 < a2 with u = h on the boundary. The Laplace equation in polar coordinates (r, θ) is 1 (rur )r + (uθ )θ = 0 r A “separation of variables” trial solution u(r, θ) = R(r)Θ(θ) gives 1 1 R00 Θ + R0 Θ + 2 RΘ00 = 0 r r which can be rearranged to Θ00 r2 R00 rR0 + =− . R R Θ Equating both sides to the separation constant λ, we are left with two ODEs, Θ00 + λΘ = 0

(1)

and r2 R00 + rR0 − λR = 0.

(2) √ For λ 6= 0, the general solution of (1) is Θ(θ) = A cos λθ + B sin λθ. Substituting this into the periodic boundary conditions Θ(0) = Θ(2π) and Θ0 (0) = Θ0 (2π) gives the homogeneous equations √ √ A(−1 + cos(2π λ)) + B sin(2π λ) = 0 √ √ √ √ −A λ sin(2π λ) + B λ(−1 + cos(2π λ)) = 0 √

This system of equations has a nontrivial solution (that is, a solution other than A =√B = 0) if the√determinant of the coefficient matrix is zero. The determinant is 2 λ(1−cos(2π λ)), and the nonzero values of λ that give a nontrivial solution are λn = n2 with n = 1, 2, . . . For λ = 0, the general solution of (1) is Θ(θ) = A + Bθ, and the only nontrivial periodic solution is Θ(θ) = A (nonzero constant). Consider now the equation (2). For λ 6= 0, a trial solution of the form R(r) = α r gives (α2 − λ)rα = 0, √ which implies the solutions α = ± λ = ±n. For λ = 0, (2) reduces to r(rR0 )0 = 0, which has the general solution R(r) = c1 ln r + c2 .

1

Writing the solution as a linear combination of the solutions found above, we have X u(r, θ) = c1 ln r + c2 + (An cos nθ + Bn sin nθ)rn n≥1

+

X

˜n sin nθ)r−n (A˜n cos nθ + B

(3)

n≥1

˜n = 0, To ensure continuity and boundedness at the origin, we set c1 = A˜n = B leaving us with X u(r, θ) = 12 A0 + (An cos nθ + Bn sin nθ)rn n≥1

where the constant term c2 has been renamed 12 A0 . At the boundary r = a we have X (An cos nθ + Bn sin nθ)an = h(θ) u(a, θ) = 12 A0 + n≥1

and the coefficents are determined by equating An an and Bn an with the Fourier series coefficients of h:  Z π X  1 Z π 1 0 0 0 0 0 u(r, θ) = h(θ ) dθ + h(θ ) cos nθ dθ cos nθ 2π −π π −π n≥1     Z π r n 1 0 0 0 h(θ ) sin nθ dθ sin nθ (4) + π −π a Example 1 Find the steady-state temperature in a long cylinder of radius a if the upper half is kept at u = 100 and the lower half is kept at u = 0. Solution. The boundary function is h(θ) = 50 + 50f (θ) where 1   −1 −π < θ < 0 θ 0 θ=0 f (θ) = -π π  1 0 < θ < π. -1

In section 8 we derived the Fourier series X 2(1 − (−1)n ) f (θ) = sin(nθ), nπ n≥1 so the PDE solution can be written directly as u(r, θ) = 50+50

X 2(1 − (−1)n )  r n n≥1

nπ

a

sin(nθ)

Plotting the sum of a large number of terms gives the above figure. 2 The series solution (4) for the Dirichlet problem in the disk can be written as " # Z π X  r n 1 0 0 u(r, θ) = h(θ ) 1 + 2 cos n(θ − θ ) dθ0 (5) 2π −π a n≥1 2

0

Letting z = ar ei(θ−θ ) , the term in brackets can be written X z z¯ 1 − z z¯ 1 − |z|2 n n 1+ z + z¯ = 1 + + = = 1 − z 1 − z¯ (1 − z)(1 − z¯) 1 + |z|2 − (z + z¯) n≥1 =

1 − (r/a)2 1 + (r/a)2 − 2(r/a) cos(θ − θ0 )

Substituting this into (5) gives the Poisson integral formula   Z π a2 − r 2 1 0 h(θ ) 2 u(r, θ) = dθ0 . 2π −π a + r2 − 2ar cos(θ − θ0 ) One application of this formula is to provide derivation of the mean R π an0 alternative 1 0 value property: the value at r = 0 is 2π −π h(θ ) dθ , which is the average of u on the circumference r = a.

10.2

Dirichlet Problem in a rectangle

Consider the two dimensional Laplace equation ∆u = 0 in the rectangle (0, a) × (0, b) with u(x, b) = g(x) and u = 0 on the rest of the boundary. Substituting a solution of the form u(x, y) = X(x)Y (y) into uxx + uyy = 0 gives

y b u=0

X 00 Y + XY 00 = 0,

u=g

x u=0

which can be rearranged to

a

Y 00 X 00 = . −X Y Equating both sides to the separation constant λ, we are left with two ODEs, X 00 + λX = 0

(6)

and Y 00 − λY = 0.

(7) √ √ For λ 6= 0, the general solution of (6) is X = A cos λx + B sin λx. Substituting this √ into the boundary conditions X(0) = 0 and X(a) = 0 gives A = 0 and B sin( λa) = 0, and the nonzero values of λ that give a nontrivial solution are λn = (nπ/a)2 with n = 1, 2, . . .. For λ = 0, the general solution of (6) is X(x) = A + Bx, and there is no nontrivial solution satisfying the boundary conditions. The general solution of (7) with λ = (nπ/a)2 is Y (y) = A cosh(nπy/a) + B sinh(nπy/a). The boundary condition Y (0) = 0 is satisfied by setting A = 0. Writing the solution as a linear combination of the solutions found above, we have X u(x, y) = Bn sin(nπx/a) sinh(nπy/a) n≥1

At the boundary y = b, this is X u(x, b) = Bn sin(nπx/a) sinh(nπb/a) = g(x). n≥1

The coefficients are obtained by equating Bn sinh(nπb/a) with the Fourier sine series coefficients of g(x):  X 2 Z a 0 0 0 sin(nπx/a) sinh(nπy/a) g(x ) sin(nπx /a) dx (8) u(x, y) = a sinh(nπb/a) 0 n≥1 3

u=0

uxx + uyy = 0

Example 2 Find the steady-state temperature in a long prismatic tube with square a × a cross-section if the top face is kept at u = 1 and other three faces are kept at u = 0. Solution. The boundary function is g(x) = 1, whose Fourier sine series was found in section 8 to be g(x) =

X 2(1 − (−1)n )

sin(nπx/a).

nπ

n≥1

The Dirichlet problem solution is then u(x, y) =

X 2(1 − (−1)n ) sin(nπx/a) sinh(nπy/a) nπ

n≥1

sinh(nπ)

Plotting the sum of a large number of terms gives the above figure. 2 The solution for the general Dirichlet problem in the rectangle is found by superposition of (8) and solutions of similar problems: 0 0

• For u(x, 0) = h(x), and u = 0 elsewhere

0,

the solution is

h

u(x, y) =

X 2 Z n≥1

a

a 0

0

h(x ) sin(nπx /a) dx

0



0

sin(nπx/a) sinh(nπ(b − y)/a) sinh(nπb/a)

This is found by replacing y by b − y [the laplacian is invariant to reflection] and g by h in (8). 0 k,

0

• For u(a, y) = k(y), and u = 0 elsewhere

the solution is

0

u(x, y) =

X 2 Z n≥1

b

b 0

0

k(y ) sin(nπy /b) dy

0



0

sin(nπy/b) sinh(nπx/b) sinh(nπa/b)

which is found by swapping x ↔ y and a ↔ b and replacing g by k in (8). 0

• for u(0, y) = j(y), and u = 0 elsewhere

j

0,

the solution is

0

u(x, y) =

X 2 Z n≥1

10.3

b

b 0

0

j(y ) sin(nπy /b) dy

0

0



sin(nπy/b) sinh(nπ(a − x)/b) . sinh(nπa/b)

Dirichlet-Neumann Problem in a wedge

Consider the two dimensional Laplace equation in the sector {(r, θ) : 0 < θ < γ, r < a}, with boundary conditions u = 0 on the rays θ = 0 and θ = γ and a Neumann condition ur = h on the perimeter r = a. Assuming a solution of the form u(r, θ) = R(r)Θ(θ) and proceeding as for the disk, we 4

ur = h

u=0 γ u=0

obtain the ODEs (1) and (2). From the boundary conditions Θ(0) = Θ(γ) = 0 we obtain Θ(θ) = sin(nπθ/γ) with n = 1, 2, . . . and R(r) = r±nπ/γ . Writing the solution as a linear combination, we have X ˜n sin(nπθ/γ)r−nπ/γ u(r, θ) = Bn sin(nπθ/γ)rnπ/γ + B n≥1

˜n = 0, leaving To ensure boundedness at the origin, we set B X u(r, θ) = Bn sin(nπθ/γ)rnπ/γ n≥1

Substituting this into the Neumann boundary condition gives X nπ nπ h(θ) = ur (a, θ) = Bn sin(nπθ/γ)a γ −1 γ n≥1 nπ

and the coefficients are determined by equating nπ Bn a γ −1 with the Fourier sine γ coefficients of h(θ), leading finally to   r nπ/γ X γa  2 Z γ 0 0 0 u(r, θ) = h(θ ) sin(nπθ /γ) dθ sin(nπθ/γ) nπ γ 0 a n≥1 π

having γ > π, the term r γ −1 in the first term of In a nonconvex sector the series for ur is unbounded as r → 0, and so this solution does not, strictly speaking, satisfy the PDE — in fact, the problem does not have a solution (having second derivatives continuous up to the boundary) in this case. Example 3 Solve ∆u = 0 in a γ = π/2 sector with u = 0 on the rays and ur = 1 on the perimeter r = a. Solution. We have Z 2(1 − (−1)n ) 2 γ h(θ0 ) sin(nπθ0 /γ) dθ0 = γ 0 nπ and so u(r, θ) =

X (1 − (−1)n )a n≥1

n2 π

sin(2nθ)

 r 2n a

.

Notice that Gibbs oscillation is not visible in the plot of the solution surface. 2

10.4

Dirichlet Problem in the region outside a circle

Consider the two dimensional Laplace equation in the region x2 + y 2 > a2 with u = h on the boundary and u bounded at infinity. We can proceed as in section 10.1 up to formula (3). To ensure the solution is bounded at infinity, we set c1 , An and Bn to zero, leaving X ˜n sin nθ)r−n u(r, θ) = 21 A˜0 + (A˜n cos nθ + B n≥1

5

Imposing the boundary condition on the circle perimeter gives X ˜n sin nθ)a−n h(θ) = u(a, θ) = 21 A˜0 + (A˜n cos nθ + B n≥1

˜n to the Fourier The coefficients are determined by equating a−n A˜n and and a−n B coefficients of h, leading to  Z π X  1 Z π 1 0 0 0 0 0 h(θ ) dθ + h(θ ) cos nθ dθ cos nθ u(r, θ) = 2π −π π −π n≥1  Z π    1 a n 0 0 0 + h(θ ) sin nθ dθ sin nθ . (9) π −π r An alternative derivation is to use the change of variables r0 = a2 /r and u0 (r0 , θ) = u(a2 /r0 , θ). By the chain rule, u0r0 = and

∂u0 ∂u ∂r a2 = = − ur ∂r0 ∂r ∂r0 (r0 )2

 r 4 1 r ∂r 1 0 0 1 a2 (r ur0 )r0 = 0 (− 0 ur )r0 = 2 (−rur )r 0 = (rur )r r0 r r a ∂r a r

and so ∆0 u0 =

 r 4 1 0 0 1 0 0 + (r u (u ∆u. ) ) = 0 r θ r r0 (r0 )2 θ a

Thus, if ∆0 u0 = 0 inside the circle then ∆u = 0 outside it. The solution inside the circle is (4) written with u0 and r0 , and applying the change of variables to this solution gives (9). Similarly, the Poisson integral formula for the disk interior is transformed to   Z π r 2 − a2 1 0 h(θ ) 2 dθ0 . u(r, θ) = 2π −π a + r2 − 2ar cos(θ − θ0 ) for the exterior. Example 4: Stationary flow past a circular cylinder The steady-state twodimensional velocity field for the irrotational flow of an incompressible inviscid constant-density fluid is given by ψy i−ψx j, where ψ is a harmonic function called the stream function. Because the velocity is orthogonal to ∇ψ, the lines of constant ψ (called streamlines) are tangential to the velocity field. Find the streamlines for flow past a long circular cylinder of radius a whose axis is the z axis, assuming the flow far from the cylinder to be constant in the x direction, that is, ψ = U y. Solution. The stream function satisfies the two dimensional Laplace equation on the exterior of the circle r = a. The radial component of the velocity is zero on the cylinder boundary, that is, the cylinder boundary is a streamline, which gives the boundary condition ψ =constant (say, zero) for r = a. Substituting this into the general solution (3) we find c1 ln a + c2 = 0 An an + A˜n a−n = 0 ˜n a−n = 0 Bn an + B 6

so that   r X n a2n ψ(r, θ) = c1 ln + r − n (An cos nθ + Bn sin nθ) a n≥1 r In order to satisfy the condition ψ = U y = U r sin θ at r → ∞ we set B1 = U and the remaining An and Bn coefficients to zero, leaving   a2 r ψ(r, θ) = c1 ln + U 1 − 2 r sin θ. a r Here are streamlines for various c1 values, which correspond to different cylinder clockwise rotation speeds:

c =0

c = Ua

c1 = 2Ua

c1 = 2.1Ua

1

1

The closely spaced streamlines correspond to regions of low pressure and indicate the presence of a net “lift” force in the y direction (Magnus effect).

7