Today’s Problem

Using only Python, solve the 2D Laplace equation

on a semi-infinite rectangular domain

with the boundary conditions

and

as well as

Solution with Python

The solution uses the separation of variables method which can be applied to a wide range of linear partial differential equations. That is, we make the ansatz

where X and Y are functions of one variable, only. Note that this ansatz will almost never give the solution directly! Instead, it will give us a base of solutions and as the Laplace equation is linear, we form linear combinations of our base solutions to match the boundary conditions. The linear combination, of course, is no longer a product of the form (1).

First, we define our symbols

from sympy import *

x, y = symbols("x, y")
u = Function("u")
X, Y = symbols("X, Y", cls=Function)

We also define a convenience function apply which allows to apply a given operation to both sides of an equation:

def apply(eq, func, *args):
    return Eq(
        func(eq.lhs, *args), 
        func(eq.rhs, *args)
    )

Now we formulate the Laplace equation:

EQ1 = Eq(
    diff(u(x, y), (x, 2)) + diff(u(x, y), (y, 2)),
    0
)

Substitute our separation ansatz (1):

EQ1.subs(u(x, y), X(x) * Y(y))

Note that the derivative with respect to y is standing in front of the whole product X(x)Y(y). This means that the derivative has not yet been executed. We can do it by calling, well, doit :

_.doit()

The derivative operator has been pulled through the factors that were not affected by the derivative. And since the remaining derivative only acts on a function of a single variable, sympy has replaced the partial derivative with a total derivative.

The next step in the separation method is to divide by the product X(x)Y(y):

apply(_, Mul, 1 / (X(x) * Y(y)))

We can separate that into two fractions by calling expand which works directly on equation objects so that we don’t need our custom function apply :

expand(_)

We put the two terms on different sides of the equation adding the negative of the second ( args[1] ) term of the left-hand side:

apply(_, Add, - _.lhs.args[1])

We will save that equation for later reference:

EQ2 = _

Now comes the crucial point of the separation method. The left-hand side of the equation only depends on x, the right-hand side only on y. Select a y and keep it fixed. Then the right side is just a number, say − k ². x is still allowed to vary. So for any x, the left-hand side is equal to a fixed number:

k = symbols("k")

Eq(EQ2.lhs, -k**2)

apply(_, Mul, X(x))

EQ3 = _

This is just the ordinary differential equation of a harmonic oscillator! We can immediately write down the solution for X(x):

var("A, B")

Eq(X(x), A * sin(k*x) + B * cos(k*x))

Instead of using symbols or Symbol , we have used var to inject the new symbols directly into the global namespace of our notebook.

The boundary conditions at x=0 requires that B=0:

_.subs(B, 0)

EQ4 = _

And the boundary condition at x=a gives a constraint on k:

var("a")

Eq(sin(k*a), 0)

solve(_, k*a)

Out: [0, pi]

This means that ka=0 and ka=π are allowed solutions for ka, but of course, also solutions shifted by integer multiples of 2π. So we have

var("n", integer=True)

Eq(k*a, n * pi)

solve(_, k)

Out: [pi*n/a]

Eq(k, _[0])

EQ_k = _

So we can eliminate the real k from our expressions and replace it, effectively, with the integer n:

EQ4.subs(k, _.rhs)

EQ5 = _

But is n=0 really a solution? No, because this would imply u(x,y)=0 and that would violate the boundary condition u(x,0)=u0.

Now let’s turn our attention to the solution of Y(y). We have

Eq(EQ2.rhs, -k**2)

apply(_, Mul, -Y(y))

This is the hyperbolic counterpart of the harmonic oscillator, so we can immediately write down, but this time in the exponential form:

var("C, D")

Eq(Y(y), C * exp(k*y) + D * exp(-k*y))

As the boundary condition for y→∞ requires the solution to go to zero, we can discard the term with exp⁡(ky):

_.subs({C: 0})

EQ6 = _

Now we can build our preliminary solution (one boundary condition is yet to be evaluated!). Note that the constants A and D above may depend on k and therefore also on n. So we attach indices. Then our preliminary solution is

A_n, D_n = symbols("A_n, D_n")

Eq(
    u(x, y),
    Sum(
        EQ5.rhs * EQ6.rhs,
        (n, 1, oo)
    ).subs({A: A_n, D: D_n, k: EQ_k.rhs})
)

But we don’t need the constants An and Dn separately. So we can introduce constants bn=An Dn:

b = IndexedBase("b")

_.subs(A_n*D_n, b[n])

EQ7 = _

The only boundary condition that remains to be evaluated is the one for y=0. This will give us the remaining coefficients bn.

u0 = Symbol("u_0")

Eq(u0, EQ7.rhs.subs(y, 0))

How can we solve this equation for bn? That’s easy if you know that the sines constitute a set of orthogonal functions. That is, their scalar product satisfies

So we simply multiply the equation by a suitable sine:

var("m", integer=True, positive=True)

apply(_, Mul, sin(pi*m*x/a))

and integrate over x:

apply(_, Integral, (x, 0, a))

EQ8 = _

On the right-hand side we can swap integration and summation:

Eq(EQ8.lhs, 
    Sum(
        b[n] * 
        Integral(sin(pi*m*x/a) * sin(pi*n*x/a), (x, 0, a)),
        (n, 1, oo)
    )
)

But due to the orthogonality property (2) of the sine function, all of the terms in the sum vanish except the one for n=m. So we have

Eq(_.lhs, _.rhs.args[0].subs(n, m))

EQ9 = _

Evaluate the integrals:

EQ9.doit()

But we know that m is positive so we can use this as an assumption to simplify with refine :

apply(_, refine,  Q.nonzero(pi*m/a))

solve(_, b[m])[0]

Eq(b[n], _.subs(m, n))

This completes our solution:

EQ7.replace(b[n], _.rhs)

SOL = _

This solution is exact, but of course, due to the infinite sum, we cannot evaluate it exactly for given x and y. But we can get arbitrarily good approximations by cutting off the sum at some point. For example, with 30 terms we get:

Sum(SOL.rhs.args[0], (n, 1, 30))

Take for instance u0=100 and a=10:

SOL_APPROX = _.subs({u0: 100, a: 10})

Then we can easily plot that approximate solution using plot_contour from sympy plotting backends spb :

from spb import *
plot_contour(SOL_APPROX, (x, 0, 10), (y, 0, 8));

If you need a better approximation, just take more terms.

This concludes the solution for the Laplace equation.