Laplace Equation in 2D on a Square Grid

The Laplace equation $ \nabla^2 u = 0 $ in rectilinear coordinates in 2D is

$$\frac{\partial^2}{\partial x^2} u(x,y) + \frac{\partial^2}{\partial y^2} u(x,y) = 0$$

We will solve it on a square region $ [0,1]\times[0,1] $ with boundary conditions $ u(0,y) = u(1,y) = 0 $ and $ u(x,0) = u(x,1) = F(x) $, where

$$ F(x) = f_1(x) - \frac 1 9 f_3(x) + \frac 1 {25} f_5(x) \quad \text{and} \quad f_n(x) = \sin(n \pi x)$$

Using separation of variables, let $ u(x,y) = f(x) \, g(y) $. Then

$$ f''(x) = \lambda f(x) \quad \text{and} \quad g''(y) = -\lambda g(y) $$

We find $ \lambda_n = - (n\pi)^2 $, the general solution to $ \nabla^2 u = 0 $ with $ u(0,y) = u(1,y) = 0 $ is

$$ u(x,y) = \sum_{n=1}^{\infty} f_n(x)\left(a_n e^{n\pi y} + b_n e^{-n\pi y}\right) $$

The exponential factors describe how each sine mode extends into the interior of the square.

Applying boundary conditions $ u(x,0) = u(x,1) = F(x) $ gives

$$ u(x,0) = \sum_{n=1}^{\infty} (a_n + b_n) f_n(x) = F(x) $$ $$ u(x,1) = \sum_{n=1}^{\infty} \left(a_n e^{n\pi} + b_n e^{-n\pi}\right) f_n(x) = F(x) $$

The coefficients of the sine modes must now be matched. For example, for $ n = 3 $, $ F(x) $ contains a term $ -\frac{1}{9}f_3. $ Therefore we must have

$$ a_3 + b_3 = -\frac{1}{9} $$ $$ a_3 e^{3\pi} + b_3 e^{-3\pi} = -\frac{1}{9} $$

This can be solved for $ a_3 $ and $ b_3 $. We get

$$ a_3 = -\frac{1}{9} \frac{e^{3\pi} - 1}{e^{9\pi} - 1} = -\frac{1}{9} \frac{1}{e^{3\pi} + 1} $$

and

$$ b_3 = -\frac{1}{9} \frac{e^{-3\pi} - 1}{e^{-9\pi} - 1} = -\frac{1}{9} \frac{1}{e^{-3\pi} + 1} $$ $$ = -\frac{1}{9} \frac{e^{3\pi}}{e^{3\pi} + 1} $$

Similarly, we can compute the other non-zero coefficients. A shortcut involves realizing that the boundary conditions are symmetric under reflecting $ y \to 1-y $. A symmetric combination of exponential functions is called the hyperbolic cosine:

$$ \cosh z = \frac{1}{2}\left(e^z + e^{-z}\right) $$

The function $ \cosh\left(y-\frac{1}{2}\right) $ is shifted so that it is clearly symmetric under $ y \to 1-y $, since

$$ \cosh\left(1-y-\tfrac{1}{2}\right) = \cosh\left(\tfrac{1}{2}-y\right) = \cosh\left(y-\tfrac{1}{2}\right) $$

An antisymmetric combination of exponential functions is called the hyperbolic sine:

$$ \sinh z = \frac{1}{2}\left(e^z - e^{-z}\right) $$

The function $ \sinh y $ is antisymmetric and equals zero at $ y = 0 $. The function $ \sinh\left(y-\frac{1}{2}\right) $ is zero at $ y = \frac{1}{2} $ and is antisymmetric under $ y \to 1-y $, since

$$ \sinh\left(1-y-\tfrac{1}{2}\right) = \sinh\left(\tfrac{1}{2}-y\right) = -\sinh\left(y-\tfrac{1}{2}\right) $$

Using these functions, we can write a different formula for the general solution.

$$ u(x,y) = \sum_{n=1}^{\infty} f_n(x)\left(c_n \cosh\big(n\pi(y - \tfrac{1}{2})\big) + d_n \sinh\big(n\pi(y - \tfrac{1}{2})\big)\right) $$

Imposing symmetry under $ y \to 1 - y $ gives $ d_n = 0 $ (only the symmetric part survives because the boundary conditions themselves are symmetric). Impose boundary conditions

$$ u(x,0) = u(x,1) = \sum_{n=1}^{\infty} c_n \cosh \left( \frac{n\pi}{2} \right) f_n(x) = F(x) $$

to read off that $ c_n $ is the coefficient of $ f_n $ in $ F(x) $ divided by $ \cosh ( \tfrac{n\pi}{2} ) $. This shortcut involves a lot less algebra, but more deep thinking.

The derivation above was carried out for the square domain $ [0,1]\times[0,1] $. In the jupyter notebook environment linked below, we extend the same idea to a rectangular domain $ [0,1]\times[0,L] $. The only structural change is that the vertical midpoint moves from $ y=\tfrac{1}{2} $ to $ y=\tfrac{L}{2} $, so the symmetric hyperbolic cosine factor becomes

$$ \cosh\left(n\pi \left(y-\tfrac{L}{2}\right)\right) \quad \text{normalized by} \quad \cosh\left(\tfrac{n\pi L}{2}\right) $$

This normalization ensures that each sine mode has the correct amplitude at both $ y=0 $ and $ y=L $. By plotting different values of $L$, we can see how the height of the rectangle changes the way the boundary data extends into the interior.

The plots below show the resulting solution for three different rectangle heights. When $ L $ is small, the two horizontal boundaries are close together, so the boundary data strongly influences the whole domain. As $ L $ increases, the boundary values become more localized near $ y=0 $ and $ y=L $, with a flatter region forming near the middle.

To see the code for these plots, click the link below:

Open Notebook Environment