Introduction Notes Theory HOWTO Examples Engineering Error Questions Matlab Maple

In this topic, we look at linear elliptic partial-differential equations (PDEs) and examine how we can solve the when subject to Dirichlet boundary conditions.

Theory

Recall that u_x(x, y) is a convenient short-hand notation to represent the first partial derivative of u(x, y) with respect to x. Given the general linear 2nd-order partial-differential equation in two variables:

a(x, y) u_{x, x}(x, y) + b(x, y) u_{x, y}(x, y) + c(x, y) u_{y, y}(x, y) +
d(x, y) u_x(x, y) + e(x, y) u_y(x, y) +
f(x, y) u(x, y) = g(x, y)

Such a PDE is termed elliptical if a(x, y) c(x, y) − b(x, y)². In the case of constant coefficients, this simplifies to ac − b². As with the finite-difference method, we can replace each of the partial derivatives with their centred divided-difference formulae, however, we will focus on two forms of this equation, namely, Poisson's equation:

u_{x, x}(x, y) + u_{y, y}(x, y) = g

and the special case when g = 0, Laplace's equation:

u_{x, x}(x, y) + u_{y, y}(x, y) = 0

We will solve Laplace's equation over a connected region R given the value of u(x, y) on the boundary ∂R. Such conditions are referred to as Dirichlet boundary conditions. We will denote such a boundary-value function by ∂u(x, y).

For simplicity, we will start by considering a rectangular region R = [x_a, x_b] × [y_a, y_b], as is shown in Figure 1 with the boundary marked in red.

Figure 1. The square region R = [x_a, x_b] × [y_a, y_b].

We will assume that x_b − x_a = nh and y_b − y_a = mh for appropriate integers n and m, and therefore, we can divide the intervals [x_a, x_b] and [y_a, y_b] into n and m sub-intervals, respectively, by defining:

x_i = x_a + ih
y_i = y_a + jh

for i = 0, 1, 2, ..., n and j = 0, 1, 2, ..., m. We will try to approximate the value of u(x, y) at each of the points (x_i, y_i), and thus, we represent these approximation by u_i, j ≈ u(x_i, y_i).

These points and the corresponding approximate values are given in Figure 2. The boundary values are once again marked in red.

Figure 2. The subdivided region. Click to enlarge.

This region has four boundaries, and from the Dirichlet conditions, we have

u_{i, 0} = ∂u(x_i, y_a)
u_{i, m} = ∂u(x_i, y_b)
u_{0, j} = ∂u(x_a, y_j)
u_{n, j} = ∂u(x_b,y_j)

Therefore, the values at each of the boundary points (red) shown in Figure 2 are known.

The Finite-Difference Equation

The next step is to convert the partial-differential equation into a finite-difference equation. From Topic 13, we can replace each second derivative with its centred divided-difference approximation.

Multiplying by h², we get:

If g = 0 (Laplace's equation), we note that this simplified to

which is independent of scaling.

Developing the System of Equations

Now, for each interior point (x_i, y_j) where i = 1, ..., n − 1 and j = 1, ..., m − 1, we can write down the finite difference equation corresponding to those points. For example, given i = 1 and j = 1, we have:

u_{2, 1} + u_{0, 1} + u_{1, 2} + u_{1, 0} − 4 u_{1, 1} = gh²

In this case, two of the points are on the border (u_{0, 1} and u_{1, 0}). From our boundary condition, we can find the value at these point:

u_{2, 1} + u_{1, 2} − 4 u_{1, 1} = gh² − ∂u(x₀, y₁) − ∂u( x₁, x₀)

This is shown graphically in Figure 3, the two border points (in red) were moved to the right-hand side of the equation.

Figure 3. Evaluating the divided-difference formula at (i, j) = (1, 1).

If we repeat this process at i = 2 and j = 1, we have:

u_{3, 1} + u_{1, 1} + u_{2, 2} + u_{2, 0} − 4 u_{2, 1} = gh²

In this case, only one of the points is on the border (u_{2, 0}). From our boundary condition, we can find the value at this point:

u_{3, 1} + u_{1, 1} + u_{2, 2} − 4 u_{2, 1} = gh² − ∂u(x₂, y₀)

This is demonstrated in Figure 5.

Figure 4. Evaluating the divided-difference formula at (i, j) = (2, 1).

Repeating this process at each interior point, we come up with a system of (n − 1)(m − 1) linear equations in (n − 1)(m − 1) unknowns. Having such a system, we can now solve for the interior values.

Approximating Intermediate Values using Interpolation

We can now use interpolation to estimate intermediate values. For example, suppose that a point (x, y) ∈ [x_{i − 1}, x_i] × [x_{i − 1}, x_i], then we can find the interpolating bivariate polynomial of the form c₁xy + c₂x + c₃y + c₄ by solving

and then evaluating the polynomial at the point (x, y).

HOWTO

Problem

Approximate the solution of a Poisson equation

u_xx(x, y) + u_xx(x, y) = g

on a rectangular region R = [x_a, x_b] × [y_a, y_b] with Dirichlet boundary conditions given by a function ∂u(x, y).

Assumptions

The function ∂u(x, y) must be piecewise continuous. We will also assume that x_b − x_a = hn and y_b − y_a = hm for some integer values of n and m.

Tools

We will use the centred divided-difference formula for the partial derivatives and linear algebra.

Process

Divide the interval [x_a, x_b] into n sub-intervals by setting x_i = x_a + ih for i = 0, 1, 2, ..., n and y_i = y_a + jh for j = 0, 1, 2, ..., m. Let u_{i, j} represent the approximation of the solution u(x_i, y_j).

At each interior point (x_i, y_j) where i = 1, 2, ..., n − 1 and j = 1, 2, ..., m − 1, evaluate the finite-difference equation:

u_{i + 1, j} + u_{i − 1, j} + u_{i, j + 1} + u_{i, j − 1} − 4 u_{i, j} = gh²

This defines a system of (n − 1)(m − 1) linear equations and (n − 1)(m − 1) unknowns. We can now solve these for the approximations u_i, j.

Examples

Example 1

Suppose we are solving Laplace's equation on [0, 1] × [0, 1] with the boundary condition defined by

using h = 0.25.

Ordering the interior points by (1, 1), (1, 2), (1, 3), (2, 1), (2, 2), (2, 3), (3, 1), (3, 2), (3, 3), we associate the following points. Thus, we can associate each of the columns and rows of the matrix as is shown in the matrix<

If we evaluate Laplace's equation at the first point, (1, 1), we get:

u_{2, 1} + u_{0, 1} + u_{1, 2} + u_{1, 0} − 4 u_{1, 1} = 0

Two of these points are on the boundary, so we can evaluate the boundary function ∂u at these to get ∂u(0.25, 0) = ∂u(0, 0.25) = 0, and thus, our linear equation simplifies to:

u_{2, 1} + u_{1, 2} − 4 u_{1, 1} = 0

Thus, we fill in the first row of the matirx.

Choosing another point, (3, 2), we get:

u_{4, 2} + u_{2, 2} + u_{3, 3} + u_{3, 1} − 4 u_{3, 2} = 0

This has only one neighbouring boundary-value point where ∂u(1, 0.5) = 1, and therefore, the equation simplifies to:

u_{2, 2} + u_{3, 3} + u_{3, 1} − 4 u_{3, 2} = -1

This defines the 8th row of the matrix, which, if we continue to fill it in, we would get:

The red entries in the column vector indicate those entries which may have been effected by ∂u on the boundary.

Solving this system of equations yields the vector u = (1/7, 2/7, 1/2, 2/7, 1/2, 5/7, 1/2, 5/7, 6/7)^T. If we plot these points, we get:

Figure 1. The approximation of the solution to Laplace's equation shown in Example 1 with h = 0.25.

If we try this again with h = 0.05, we get the approximations shown in Figure 2.

Figure 2. The approximation of the solution to Laplace's equation shown in Example 1 with h = 0.05.

If we wanted a better approximation, we could use a smaller value of h.

Figure 3. The approximation of the solution to Laplace's equation shown in Example 1 with h = 0.01.

Example 2

Explore what happens when we solve Poisson's equation

for the same boundary conditions as given in Example 1 for values of g between 0 and 4.

u_xx + u_yy = g

Using n = m = 32, Figure 4 shows the approximations for values of g starting with Laplace's equation and going to g = 4.

Figure 4. The solutions to the Poisson equation for values of g ∈ [0, 4].

What Poisson's equation is dictating is that locally, the solution will look like x² + y².

Example 3

Let the boundary condition in Example 1 be replaced by the function cos(π(x + y)). Solve Laplace's equation with this boundary condition.

Using the same ordering of the interior points and using the approximation cos(π/4) ≈ 0.7071, we get the system of linear equations

Solving this system of linear equations, we get the vector

u = 0.35355339, 0, -0.35355339, 0, 0, 0, -0.35355339, 0, 0.35355339)^T

Plotting these points yields the image in Figure 5.

Figure 5. Solving Laplace's equation with sinusoidal boundary conditions.

Using a larger number of points (h = 0.01), we get the approximation shown in Figure 6. Note that the solution looks a lot like the bubble which would result if the wire frame (the boundary) was dipped into soap water.

Figure 6. Solving Laplace's equation with sinusoidal boundary conditions with h = 0.01.

Engineering

Under static conditions, electric fields obey Poisson's equation. On the large scale in a vacuum, the fields may be approximated by solutions to Laplace's equation.

Error

To be completed later.

Questions

Question 1

Approximate the solution to Laplace's equation on [0, 1] × [0, 1] with the Dirichlet boundary condition:

with h = 0.25.

Answer:

u(0.25, 0.25) ≈ 0.5
u(0.25, 0.5) ≈ 0.625
u(0.25, 0.75) ≈ 0.5
u(0.5, 0.25) ≈ 0.375
u(0.5, 0.5) ≈ 0.5
u(0.5, 0.75) ≈ 0.375
u(0.75, 0.25) ≈ 0.5
u(0.75, 0.5) ≈ 0.625
u(0.75, 0.75) ≈ 0.5

Question 2

Given the same boundary conditions as in Question 1, approximate the solution to Poisson's equation with g = 16.

Answer:

u(0.25, 0.25) ≈ -0.1875
u(0.25, 0.5) ≈ -0.25
u(0.25, 0.75) ≈ -0.1875
u(0.5, 0.25) ≈ -0.5
u(0.5, 0.5) ≈ -0.625
u(0.5, 0.75) ≈ -0.5
u(0.75, 0.25) ≈ -0.1875
u(0.75, 0.5) ≈ -0.25
u(0.75, 0.75) ≈ -0.1875

Question 3

Approximate the solution to Poisson's equation u_xx + u_yy = 1 given the boundary condition ∂u = 0 on the square [3, 5] × [4, 6] with h = 0.25.

Hint: use symmetry to reduce this to a system with 10 equations. Be sure to use Matlab.

Answer:

u(3.25, 4.25) ≈ -0.017779
u(3.25, 4.5) ≈ -0.027746
u(3.25, 4.75) ≈ -0.032916
u(3.25, 5) ≈ -0.034524
u(3.5, 4.5) ≈ -0.044663
u(3.5, 4.75) ≈ -0.053768
u(3.5, 5) ≈ -0.056641
u(3.75, 4.75) ≈ -0.065229
u(3.75, 5) ≈ -0.068876
u(4, 5) ≈ -0.072783

Question 4

Given the boundary conditions in Question 3, what is the solution if we were solving Laplace's equation?

Matlab

Maple

Solving Laplace's equation on [0, 1] × [0, 1] where the boundary value is defined by a function f(x, y) may be done as follows:

N :=32:                          # the number of intervals
f := (x, y) -> sin(4.0*x*y):

# assign the boundary values on [0, 1] x [0, 1]
for i from 0 to N do
    du[i,0] := f(i/N, 0);
    du[0,i] := f(0, i/N);
    du[i,N] := f(i/N, 1);
    du[N,i] := f(1, i/N);
end do:

M := Matrix( (N - 1)*(N - 1), (N - 1)*(N - 1), datatype=float[8] ):
b := Vector( (N - 1)*(N - 1), datatype=float[8] ):

for i from 1 to N - 1 do
    for j from 1 to N - 1 do
        posn := (N - 1)*(i - 1) + j;
        M[posn, posn] := -4;

        for k from -1 to 1 by 2 do
            if assigned( du[i + k, j] ) then
                b[posn] := b[posn] - du[i + k, j];
            else
                M[posn, posn + (N - 1)*k] := 1;
            end if;
        end do;

        for k from -1 to 1 by 2 do
            if assigned( du[i, j + k] ) then
                b[posn] := b[posn] - du[i, j + k];
            else
                M[posn, posn + k] := 1;
            end if;
        end do;
    end do;
end do:

u := LinearAlgebra:-LinearSolve( M, b ):   # solve the system of linear equations

# plot the result
plots[display](
    plots[pointplot3d](
        [seq( seq( [i/N, j/N, u[(i - 1)*(N - 1) + j]], i = 1..N - 1 ), j = 1..N - 1 )],
        symbol = circle, color = black
    ),
    # the boundary values
    plots[spacecurve]( [0, t, f(0, t)], t = 0..1, thickness = 2, colour = red ),
    plots[spacecurve]( [t, 0, f(t, 0)], t = 0..1, thickness = 2, colour = red ),
    plots[spacecurve]( [1, t, f(1, t)], t = 0..1, thickness = 2, colour = red ),
    plots[spacecurve]( [t, 1, f(t, 1)], t = 0..1, thickness = 2, colour = red ),
    axes = framed
);

Topic 15.2: Elliptic Partial-Differential Equations

Theory

The Finite-Difference Equation

Developing the System of Equations

Approximating Intermediate Values using Interpolation

HOWTO

Problem

Assumptions

Tools

Process

Examples

Example 1

Example 2

Example 3

Engineering

Error

Questions

Question 1

Question 2

Question 3

Question 4

Matlab

Maple