Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

By the end of this section, you will be able to:

Set up Laplace's equation $\nabla^2 u = 0$ on a rectangle with Dirichlet boundary values on all four sides, and reduce the four-boundary problem to four one-boundary problems by superposition.
Derive the building-block solution $\varphi_n(x,y) = \sin\!\big(\tfrac{n\pi x}{a}\big)\,\dfrac{\sinh(n\pi y/a)}{\sinh(n\pi b/a)}$ from separation of variables, and explain why the sign of the separation constant is forced by the boundary conditions.
Compute the Fourier-sine coefficients of the non-homogeneous boundary by hand for at least one boundary function.
Recognise the exponential decay of high-frequency modes encoded in $\sinh(n\pi b/a)$ as the smoothing mechanism of Laplace's equation.
Build a working solver in Python and PyTorch and double-check it against an analytically known boundary.
Use the interactive viewers to watch a single mode, to fit a custom top edge, and to inspect the full interior of the rectangle as you change the mode count.

The Big Picture

"A harmonic function on a region is the function that has no choice but to interpolate its boundary values as smoothly as possible."

Laplace's equation $\nabla^2 u = u_{xx} + u_{yy} = 0$ looks almost embarrassingly simple. Two second derivatives, set to zero. Yet this equation governs the steady-state temperature of a metal plate, the gravitational and electrostatic potentials in empty space, the shape of a soap film stretched on a wire frame, the velocity potential of an ideal fluid, and the equilibrium of countless other physical systems. Anywhere you find a system that has settled down and stopped changing, Laplace is hiding underneath.

In a rectangle the equation becomes one of the cleanest stories in applied mathematics. The rectangle has four straight edges, the operator is symmetric under $x \to a - x$ , and the basis we will discover — products of sines and hyperbolic sines — is so well-matched to the geometry that the entire problem reduces to one Fourier series.

The intuition in one sentence

The boundary values are projected onto a sine basis along the active edge, and each projection is carried into the interior with a known exponential decay rate — bigger frequency means faster decay away from the boundary.

The Problem We Want to Solve

Fix a rectangle $R = [0, a] \times [0, b]$ . We seek a function $u(x, y)$ defined on $R$ that satisfies

\displaystyle u_{xx}(x,y) + u_{yy}(x,y) \;=\; 0 \quad\text{inside } R

together with Dirichlet boundary conditions $u = f_i$ on each of the four sides. To keep the derivation clean we start with the simplest non-trivial case: three sides clamped to zero, one side given:

Left wall: $u(0, y) = 0$ for $0 \le y \le b$ .
Right wall: $u(a, y) = 0$ for $0 \le y \le b$ .
Bottom edge: $u(x, 0) = 0$ for $0 \le x \le a$ .
Top edge: $u(x, b) = f(x)$ — the only non-zero edge.

Solving this single-active-edge case will turn out to be enough, because of linearity: if all four edges are active, we just solve four problems of the above type and add the answers. That trick is the whole reason we focus our energy on one edge.

Why the corners are not a problem

At the corners the boundary specification can clash — $u(0, b) = 0$ from the left wall but $u(0, b) = f(0)$ from the top edge. The fix is to require $f(0) = f(a) = 0$ so the two prescriptions agree. The Fourier-sine basis we are about to use enforces this automatically — every basis function vanishes at $x = 0$ and $x = a$ .

The Intuition: Steady-State Heat in a Plate

Before any algebra, picture a thin rectangular metal plate. Three edges are held in an ice bath at $0^\circ$ . Along the top edge a heater paints a pattern $f(x)$ : hot in the middle, say, and cool near the corners. Wait long enough for the temperature inside the plate to stop changing. The resulting steady temperature $u(x, y)$ is exactly the solution we seek.

Two physical facts already tell us nearly everything about the answer:

The maximum is on the boundary. Heat flows from hot to cold; an interior hotspot would lose heat to its neighbours and cool down, contradicting steady state. The peak value of $u$ is achieved at some point on the boundary — and so is the minimum.
Information from the heater decays with depth. A sharp spike on the top edge cannot punch through to the bottom; it gets smeared out by the time it has travelled even a small fraction of the rectangle height. The finer the spike, the faster it attenuates. We will see this decay quantitatively as $e^{-n\pi y / a}$ : high frequency $n$ ⇒ shallow penetration.

Both facts will fall out of the formula. Hold this picture in your head: a heated edge whose imprint diffuses smoothly across the plate.

The Strategy: Separation of Variables

We guess that the answer is a product:

u(x, y) \;=\; X(x)\, Y(y)

This is a strong, restrictive guess: most functions of two variables are not products. The reason it works here is the equation's linearity — once we find every product solution we can add them up to build the general one. Substitute the guess into Laplace's equation:

X''(x)\, Y(y) + X(x)\, Y''(y) = 0

Divide through by $X(x) Y(y)$ (legal away from the boundary, where the product is nonzero):

\dfrac{X''(x)}{X(x)} \;=\; -\,\dfrac{Y''(y)}{Y(y)}

The left side depends only on $x$ , the right side only on $y$ . The only way these two can be equal for all $(x, y)$ is if both equal the same constant. Call that constant $-\lambda$ (the minus sign is chosen to make the next two equations look standard):

X''(x) + \lambda\, X(x) = 0, \qquad Y''(y) - \lambda\, Y(y) = 0.

The hidden art of the minus sign

Whether the constant is $+\lambda$ or $-\lambda$ looks arbitrary — but it is forced by the boundary conditions. The side walls demand $X(0) = X(a) = 0$ , a homogeneous Dirichlet problem on a finite interval, and that problem has non-trivial solutions only when $\lambda > 0$ . So we set up the equation with $X'' + \lambda X = 0$ (oscillatory in $x$ ) and the partner $Y'' - \lambda Y = 0$ (exponential in $y$ ). The geometry chooses for us.

The X-Equation: Sine Modes from the Side Walls

The X-equation $X''(x) + \lambda X(x) = 0$ with side-wall conditions $X(0) = X(a) = 0$ is the classic Sturm–Liouville eigenvalue problem. Take $\lambda = k^2$ ; the general solution is

X(x) = A \cos(k x) + B \sin(k x).

The left wall $X(0) = 0$ kills the cosine: $A = 0$ . The right wall $X(a) = B \sin(k a) = 0$ demands either $B = 0$ (trivial) or $\sin(k a) = 0$ , which forces

k = \dfrac{n\pi}{a}, \quad n = 1, 2, 3, \dots

So the allowed $X$ -functions are $\sin(n\pi x / a)$ for $n = 1, 2, 3, \dots$ , and the corresponding separation constants are $\lambda_n = (n\pi/a)^2$ . These are the only building blocks compatible with the side walls — there is no other choice.

Why we can drop B = 0

The trivial solution $X \equiv 0$ gives $u \equiv 0$ , which is the unique answer only when $f \equiv 0$ on the top edge. For any non-trivial top boundary we need at least one mode with $B \ne 0$ , so we are forced into the discrete spectrum $k = n\pi/a$ .

The Y-Equation: Hyperbolic Sines from the Bottom

With $\lambda_n = (n\pi/a)^2$ fixed by the side walls, the $Y$ -equation becomes

Y''(y) - \big(n\pi/a\big)^2 Y(y) = 0.

Now the sign of the constant matters: this is the exponential ODE, not the oscillatory one. Its general solution is

Y(y) = C\,\cosh\!\big(n\pi y/a\big) + D\,\sinh\!\big(n\pi y/a\big).

The bottom-edge condition $Y(0) = 0$ (so that $u(x, 0) = X(x) Y(0) = 0$ ) forces $C = 0$ , because $\cosh(0) = 1$ while $\sinh(0) = 0$ . We are left with

Y_n(y) = D_n\, \sinh\!\big(n\pi y/a\big).

The constant $D_n$ is absorbed into the overall coefficient of the mode (we have one free constant per $n$ ; whether we put it on $X$ or on $Y$ is bookkeeping). It is convenient to normalise the hyperbolic sine so that the mode equals $1$ at the top edge:

\widetilde{Y}_n(y) = \dfrac{\sinh(n\pi y / a)}{\sinh(n\pi b / a)}.

With this normalisation the mode profile rises from $0$ at $y = 0$ to $1$ at $y = b$ , and any leftover scaling lives in the coefficient $b_n$ we will compute from the top edge.

One Building-Block Mode

Multiplying the two pieces together, the n-th building block is

\varphi_n(x, y) = \sin\!\big(n\pi x/a\big)\,\dfrac{\sinh(n\pi y/a)}{\sinh(n\pi b/a)}.

Each $\varphi_n$ is harmonic (verify by direct substitution), satisfies all three zero boundary conditions, and equals $\sin(n\pi x/a)$ on the top edge — the n-th sine wave of the boundary basis. So our problem is reduced to matching the top edge with a sum of sines, which is exactly a Fourier sine series.

n	Top-edge shape sin(nπx/a)	Vertical decay rate from y = b
1	one bump, peak at x = a/2	slowest, e^{-π(b−y)/a}
2	one positive bump, one negative	twice as fast
3	three lobes alternating sign	three times as fast
n	n lobes with n−1 sign changes	n × baseline

Interactive: A Single Mode

Slide $n$ below. Watch two things at once: on the heatmap, the top edge gets the n-th sine wave imprinted on it; in the interior, the wiggles fade exponentially as you move downward. The little side panels show the two factors $X(x)$ and $Y(y)$ — the building blocks themselves.

What to notice while you play

For $n = 1$ , the mode penetrates almost all the way to the bottom — the slow-frequency component reaches deepest.
For $n = 6$ , the wiggles are basically invisible below $y \approx 0.4$ . High frequency, shallow imprint.
The $Y(y)$ profile is always 0 at $y = 0$ and 1 at $y = b$ — the normalization we chose so that the coefficient absorbs all the boundary scaling.

Superposition: Adding Modes to Match the Top

A single mode $\varphi_n$ can only match a top boundary that already looks like $\sin(n\pi x/a)$ . For a general $f(x)$ we need a weighted sum:

u(x, y) = \sum_{n=1}^{\infty} b_n\, \sin\!\big(n\pi x/a\big)\,\dfrac{\sinh(n\pi y/a)}{\sinh(n\pi b/a)}.

Because each $\varphi_n$ is harmonic and Laplace's equation is linear, every finite sum (and every convergent infinite sum) is also harmonic. Each term carries its share of the boundary; together they have to add up exactly to $f(x)$ on the top edge. Setting $y = b$ and using $\sinh(n\pi b/a)/\sinh(n\pi b/a) = 1$ :

f(x) = \sum_{n=1}^{\infty} b_n\, \sin\!\big(n\pi x/a\big).

That last line is the heart of the matter. The whole 2D problem reduces to a 1D Fourier sine expansion of the top-edge function. The coefficients $b_n$ we read off here immediately feed into the 2D formula above.

Reading the Top Edge as a Fourier Sine Series

How do we extract $b_n$ from $f(x)$ ? The sines $\sin(n\pi x/a)$ are orthogonal on $[0, a]$ :

\int_0^a \sin\!\big(m\pi x/a\big)\,\sin\!\big(n\pi x/a\big)\, dx = \begin{cases} 0 & m \ne n \\ a/2 & m = n. \end{cases}

Multiply the boundary identity $f(x) = \sum_n b_n \sin(n\pi x/a)$ by $\sin(m\pi x/a)$ and integrate from $0$ to $a$ . By orthogonality, every term in the sum vanishes except $n = m$ , leaving

\int_0^a f(x)\sin\!\big(m\pi x/a\big)\, dx = b_m \cdot \dfrac{a}{2}.

Solving for $b_m$ (and renaming $m \to n$ ):

b_n = \dfrac{2}{a} \int_0^a f(x)\sin\!\big(n\pi x/a\big)\, dx.

This is the only formula you need to remember. The 2D problem is completely determined by these 1D integrals along the active edge.

Interactive: Fitting the Top Boundary

Pick a top-edge function and slide $N$ . The thick black curve is the true $f(x)$ ; the red curve is the partial sum $f_N(x) = \sum_{n=1}^N b_n \sin(n\pi x)$ . The stem plot underneath shows all the coefficients $b_n$ (red = kept, grey = thrown away).

The shape of b_n tells you everything

Smooth bump: coefficients decay exponentially. Three modes essentially nail it.
Triangular pluck: coefficients decay like $1/n^2$ because the function has a corner (kink in the first derivative).
Square plateau: jumps at $x = 1/4$ and $3/4$ . Coefficients decay only like $1/n$ ; you get persistent ~9% overshoot at the jumps no matter how many modes you keep — the Gibbs phenomenon.
Single mode sin(3πx): exactly one non-zero coefficient, $b_3 = 1$ . The fit is already perfect at $N = 3$ .

The Full Solution Formula

Putting it all together:

\boxed{\; u(x,y) = \sum_{n=1}^{\infty} b_n\, \sin\!\big(n\pi x/a\big)\,\dfrac{\sinh(n\pi y/a)}{\sinh(n\pi b/a)}, \quad b_n = \dfrac{2}{a}\int_0^a f(x)\,\sin\!\big(n\pi x/a\big)\, dx. \;}

Every symbol has a meaning we have earned:

Symbol	Role	Where it came from
sin(nπx/a)	horizontal shape of mode n	side-wall boundary conditions
sinh(nπy/a)	vertical growth from bottom	bottom-edge condition Y(0) = 0
sinh(nπb/a)	normalisation at top	convenience: makes the mode equal 1 on the top
b_n	weight of mode n	Fourier inner product with f on the top edge

Interactive: The Full Rectangle Solver

Choose a top-boundary function. Slide $N$ (number of modes summed) and $b$ (rectangle height). The interior of the rectangle is the harmonic function determined by the chosen boundary. Red is positive, blue is negative, white is zero.

Three experiments to try

Pick Square plateau, set $N = 40$ . Notice the boundary fits poorly near the jumps but the interior looks perfectly smooth. Harmonic functions are instantly analytic in the interior — only the boundary values inherit the roughness.
Pick Triangular pluck, slide $b$ from $0.3$ (short rectangle) to $1.5$ (tall rectangle). The taller the rectangle, the cleaner the bottom-half interior becomes — the high-frequency boundary information cannot reach down that far.
Pick sin(3πx) at $N = 1$ . The fit is wrong (only the first mode contributes). Increase $N$ to $3$ and watch the answer snap into place — this boundary is a single building block, so three modes is overkill but two is short.

Worked Numerical Example

We commit to specific numbers and turn the crank by hand. Set up the problem so the answer is something we can recognise.

Hand-traced example: top edge

f(x) = \sin(2\pi x) + 2\sin(5\pi x)

, rectangle

a = 1, b = 1

Step 1 — Compute the coefficients. The top-edge function is already a linear combination of basis sines, so the integrals are pure orthogonality bookkeeping. Using $\int_0^1 \sin(m\pi x)\sin(n\pi x)\,dx = \tfrac{1}{2}\delta_{mn}$ :

b_n = 2\int_0^1 \big[\sin(2\pi x) + 2\sin(5\pi x)\big]\sin(n\pi x)\, dx.

For $n = 2$ only the first term contributes $(2)(1/2) = 1$ : $b_2 = 1$ . For $n = 5$ only the second term contributes $(2)(2)(1/2) = 2$ : $b_5 = 2$ . All other $b_n = 0$ .

Step 2 — Plug into the master formula with $a = b = 1$ :

u(x, y) = \sin(2\pi x)\,\dfrac{\sinh(2\pi y)}{\sinh(2\pi)} + 2\sin(5\pi x)\,\dfrac{\sinh(5\pi y)}{\sinh(5\pi)}.

Step 3 — Sanity-check the boundary.

At $y = 1$ : $\sinh(n\pi)/\sinh(n\pi) = 1$ , so $u(x, 1) = \sin(2\pi x) + 2\sin(5\pi x) = f(x)$ . ✓
At $y = 0$ : $\sinh(0) = 0$ , so $u(x, 0) = 0$ . ✓
At $x = 0$ and $x = 1$ : $\sin(2\pi \cdot 0) = \sin(5\pi \cdot 0) = 0$ , so $u(0, y) = u(1, y) = 0$ . ✓

Step 4 — Evaluate at a single point. Take $(x, y) = (1/2,\, 1/2)$ . Then $\sin(\pi) = 0$ (kills the first term) and $\sin(5\pi/2) = 1$ (full amplitude for the second).

u(1/2, 1/2) = 0 + 2 \cdot 1 \cdot \dfrac{\sinh(5\pi/2)}{\sinh(5\pi)}.

Numerically $\sinh(5\pi/2) \approx 1287.99$ and $\sinh(5\pi) \approx 3.317\times 10^{6}$ , so

u(1/2, 1/2) \approx 2 \cdot \dfrac{1287.99}{3.317\times 10^{6}} \approx 7.76 \times 10^{-4}.

Step 5 — Interpret. The top edge is order $1$ at $x = 1/2$ (the second term peaks there with value $2$ ), but by the midline $y = 1/2$ the contribution from the $n = 5$ mode has been attenuated by a factor of $\sinh(5\pi/2)/\sinh(5\pi) \approx 4 \times 10^{-4}$ . That is the price of frequency: each unit of $n$ roughly halves the depth of penetration in this rectangle.

Sanity check this with the interactive solver

Open the rectangle solver above, select sin(3πx), push $N$ to 1, change b to 1.0. The colour scale shows the same attenuation effect we just computed by hand for a different boundary — the principle is identical.

Python: Computing the Coefficients by Hand

For boundaries that are not already a sum of sines, the integral has to be done numerically. The midpoint rule with a few hundred points is enough to recover the coefficients to 4–5 digits. Here is the cleanest possible implementation, with each line annotated.

Fourier-sine coefficients of a top-edge function

🐍laplace_coefficients.py

Explanation(10)

Code(28)

1NumPy only

We are still in the 'derivation' phase, so we stay in plain NumPy. There is no neural network yet — just arrays and trig. The point of this snippet is to make the Fourier-coefficient formula b_n = (2/a) ∫_0^a f(x) sin(n π x / a) dx feel concrete by computing it with our own hands.

4Rectangle dimensions a, b

a is the width along x, b is the height along y. We chose a = 1, b = 0.8 — a slightly squat rectangle. The coefficient formula depends only on a (the width). b enters later, through the hyperbolic decay factor sinh(n π y / a) / sinh(n π b / a).

6Top-boundary function f(x)

f(x) is the value of u on the top edge y = b. Here it is a sharp Gaussian bump centred at x = 0.5 with variance 0.01. The boundary already satisfies f(0) = f(a) ≈ 0 because the bump dies off quickly — so the Fourier sine series can match it without Gibbs trouble.

EXAMPLE

f(0.5) ≈ 1.0   f(0) ≈ 3.7e-6

9Coefficient formula

This whole function is one line of calculus: b_n = (2/a) ∫_0^a f(x) sin(n π x / a) dx. The 2/a in front is the normalization of the sine basis. The integral measures how much f(x) 'looks like' the n-th sine mode. If f and sin(n π x / a) agree in shape, b_n is large.

11Midpoint rule integration

We approximate the integral with the midpoint rule: split [0, a] into M = 800 cells, evaluate f and the sine at each cell centre, multiply, sum, scale by the cell width a/M. The midpoint rule has error O(h²) for smooth integrands and is the cheapest way to get 4–5 correct digits on this kind of bump.

EXAMPLE

xs = [0.000625, 0.001875, ..., 0.999375]

12The product f(x)·sin(n π x / a)

Element-wise multiply: for each grid point we pair the value of f(x) with the value of sin(n π x / a). If both are positive at x, the product is positive; if they have opposite signs, the product is negative. Summing the products is the discrete inner product 〈f, sin_n〉.

12Why the 2 in front

The sine functions sin(n π x / a) form an orthogonal basis on [0, a] with ∫_0^a sin²(n π x / a) dx = a/2. To project f onto a unit-norm basis we divide by a/2, i.e. multiply by 2/a. That is where the 2 comes from — not magic, just normalization.

15First 6 coefficients

We compute b_1, b_2, …, b_6. Because f is symmetric about x = a/2, the even-n coefficients vanish exactly (the odd sines are symmetric, the even sines are antisymmetric — their inner product with a symmetric f is zero). You should see b_2 = b_4 = b_6 ≈ 0 when you run this.

20Reconstruct the top edge

For each x we sum b_n · sin(n π x / a) over the 6 modes. This is the partial Fourier series f_6(x). Comparing it to the true f(x) is the most honest test of 'how well is the basis tracking my boundary?' for the truncation N = 6.

23L∞ error of the partial sum

max |f(x) − f_6(x)| is the worst-case mismatch over the whole edge. For a smooth Gaussian, six modes already drop this below 1%. For a square wave with jumps, no matter how many modes you take, this number plateaus near ~9% — the Gibbs phenomenon. The decay rate of the L∞ error is what we are really paying for when we add modes.

18 lines without explanation

1import numpy as np
2
3# Rectangle [0, a] x [0, b].  Three sides clamped at u = 0.
4# Top edge carries a Gaussian bump centred at x = a/2.
5a, b = 1.0, 0.8
6
7def f(x):
8    return np.exp(-((x - 0.5) ** 2) / 0.02)
9
10# Fourier-sine coefficient:
11#   b_n = (2/a) * integral_0^a f(x) sin(n pi x / a) dx
12def b_coeff(n, M=800):
13    xs = np.linspace(0, a, M, endpoint=False) + 0.5 * a / M
14    return (2.0 / a) * np.sum(f(xs) * np.sin(n * np.pi * xs / a)) * (a / M)
15
16# Compute first 6 coefficients
17coeffs = [b_coeff(n) for n in range(1, 7)]
18for n, b_n in enumerate(coeffs, start=1):
19    print(f"b_{n} = {b_n:+.5f}")
20
21# Reconstruct the top edge from these coefficients and compare.
22xs = np.linspace(0, a, 200)
23recon = np.zeros_like(xs)
24for n, b_n in enumerate(coeffs, start=1):
25    recon += b_n * np.sin(n * np.pi * xs / a)
26
27err_inf = float(np.max(np.abs(recon - f(xs))))
28print(f"max |f(x) - reconstruction|  with 6 modes  =  {err_inf:.4f}")

Try this

Change $f(x)$ to $x(1-x)$ (a smooth parabola). Run the script. You should see only odd $n$ coefficients are non-zero (parabolic boundary is symmetric about $x = 1/2$ ) and they decay like $1/n^3$ (smoothness pays for itself).

PyTorch: A Vectorised Solver in 12 Lines

Once we have the coefficients, evaluating $u(x, y)$ on a whole grid is a single tensor operation. The same code works on CPU or GPU without modification, and it scales effortlessly to a million grid points.

Rectangle Laplace solver — one einsum, done

🐍laplace_rectangle_solver.py

Explanation(17)

Code(35)

1PyTorch import

We move to PyTorch because the whole solution can be written as a single tensor einsum once the coefficients are computed. The same machinery later plugs into a physics-informed neural network for problems where the boundary is not a clean Fourier sum.

4Rectangle dimensions

Same a = 1, b = 0.8 as the NumPy example. The variable name b in physics means 'rectangle height'; the variable name b_n means 'Fourier coefficient'. Different b's. Mathematicians have been overloading the letter b for two centuries — get used to it.

6Top-boundary function f_top(x)

We picked a boundary that is exactly a sum of two sine modes. The expected coefficients are b_3 = 1, b_5 = 0.5, and all others zero. This lets us double-check the coefficient code: if we get anything other than that, there is a bug.

9Grid sizes nx, ny

200 × 160 = 32 000 interior points. Each one needs N = 32 multiply-adds for the mode sum. In a Python loop this would be ~1 million ops and noticeable; with einsum it is a single CUDA kernel and finishes in milliseconds even on CPU.

10x grid

Equally spaced from 0 to a along the width. Shape (nx,). When we broadcast x[None, :] later, it becomes a row vector so the n index can broadcast against it cleanly.

11y grid

Equally spaced from 0 to b along the height. Shape (ny,). Note that y[0] = 0 (the bottom edge, clamped to zero) and y[-1] = b (the top edge, where u must equal f_top).

14Mode count N

32 modes is heavy overkill for this particular boundary (only two non-zero coefficients), but the code is intentionally generic — change f_top to anything else, and N = 32 will still resolve it well unless the boundary has jumps.

15Mode indices n

tensor([1.0, 2.0, ..., 32.0]). Float, not int, so that the trig and sinh calls do not promote/branch oddly. This is the n that runs through every formula in the derivation.

18Quadrature x-grid

Midpoint-rule grid for integrating along the top edge. We add 0.5 * a / 800 to shift to cell centres. 800 points × 32 modes is one (32, 800) matrix-vector op — cheap.

19Basis matrix sin_nx

sin_nx[n, j] = sin((n+1) π xq[j] / a). Shape (N, Mq). Each row is one sine basis function sampled on the quadrature grid. Computing all rows at once with broadcasting (n[:, None] vs xq[None, :]) is the whole point of switching to a tensor library.

20f sampled on quadrature grid

fq is the boundary function evaluated at every quadrature node. Shape (Mq,). We compute it once and reuse it for every coefficient — there is no point sampling f separately for each n.

21Coefficient vector b_n

Element-wise multiply fq[None, :] (1, Mq) with sin_nx (N, Mq) → broadcasts to (N, Mq). Sum along the Mq axis → (N,). Multiply by the quadrature weight a / 800 and the normalisation 2 / a. The result is the vector (b_1, b_2, …, b_N) in one line.

EXAMPLE

b_n[2] = 1.0   b_n[4] = 0.5   others ≈ 0

24Spatial sine matrix on the visualisation grid

sin_nX[n, j] = sin((n+1) π x[j] / a). Same idea as sin_nx but evaluated on the rendering grid (nx points instead of 800). Decoupling integration grid from visualisation grid is good practice — they can have very different requirements.

25Denominator sinh(n π b / a)

This is the only place where the rectangle height b enters the solution. sinh(n π b / a) grows like e^{n π b / a} / 2 for large n — meaning high-frequency modes are massively suppressed when b is large. Physically: a tall rectangle 'forgets' the top-edge wiggles long before they reach the bottom.

26Vertical decay matrix sinh_nY

sinh_nY[n, i] = sinh((n+1) π y[i] / a) / sinh((n+1) π b / a). At y = 0 this is 0 (the bottom edge condition); at y = b it is 1 (the top edge condition). In between it grows exponentially from 0 to 1 — that exponential growth is what makes the boundary energy 'reach down' into the interior.

29einsum assembly

The full solution in one tensor expression. The string 'n,nj,ni->ij' means: 'take b_n indexed by n, sin_nX indexed by (n, j), sinh_nY indexed by (n, i); multiply pointwise on n; sum over n; emit a tensor indexed by (i, j)'. That is literally the formula u(x_j, y_i) = Σ_n b_n · sin(nπx_j/a) · sinh(nπy_i/a) / sinh(nπb/a).

EXAMPLE

u.shape == (ny, nx) == (160, 200)

32Quick sanity print

We compare u[-1, nx//2] (the top edge, middle) with f_top(0.5). They should agree to a few significant digits — if not, either the coefficient code is wrong or N is too small. Always sanity-check a solver against a known boundary point.

18 lines without explanation

1import torch
2
3# Rectangle and boundary.  Top edge = sin(3 pi x) + 0.5 * sin(5 pi x).
4a, b = 1.0, 0.8
5
6def f_top(x):
7    return torch.sin(3 * torch.pi * x) + 0.5 * torch.sin(5 * torch.pi * x)
8
9# Grid: nx columns along x, ny rows along y.
10nx, ny = 200, 160
11x = torch.linspace(0.0, a, nx)            # (nx,)
12y = torch.linspace(0.0, b, ny)            # (ny,)
13
14# How many modes to sum.
15N = 32
16n = torch.arange(1, N + 1).float()        # (N,)
17
18# Fourier-sine coefficients via midpoint quadrature on the top edge.
19xq = torch.linspace(0.0, a, 800) + 0.5 * a / 800   # (Mq,)
20sin_nx = torch.sin(n[:, None] * torch.pi * xq[None, :] / a)    # (N, Mq)
21fq = f_top(xq)                                                  # (Mq,)
22b_n = (2.0 / a) * (fq[None, :] * sin_nx).sum(dim=1) * (a / 800) # (N,)
23
24# Build u(x, y) = sum_n  b_n * sin(n pi x / a) * sinh(n pi y / a) / sinh(n pi b / a)
25sin_nX = torch.sin(n[:, None] * torch.pi * x[None, :] / a)           # (N, nx)
26denom  = torch.sinh(n * torch.pi * b / a)                            # (N,)
27sinh_nY = torch.sinh(n[:, None] * torch.pi * y[None, :] / a) / denom[:, None]   # (N, ny)
28
29# Outer-product assembly:  u[i, j] = sum_n b_n * sin_nX[n, j] * sinh_nY[n, i]
30u = torch.einsum("n,nj,ni->ij", b_n, sin_nX, sinh_nY)   # (ny, nx)
31
32print("solution shape:", tuple(u.shape))
33print("max |u| in interior:", u.abs().max().item())
34print("top-edge value at x = 0.5:", u[-1, nx // 2].item(),
35      "  expected:", f_top(torch.tensor(0.5)).item())

Why einsum is the right tool here

Three rank-1 ingredients ( $b_n$ , $\sin(n\pi x/a)$ , $\sinh(n\pi y/a)/\sinh(n\pi b/a)$ ) combine to a rank-2 output via a contraction over $n$ . That is the textbook use case for $\texttt{einsum}$ . Writing it with three nested for-loops would be hundreds of times slower and ten times harder to read.

Four Active Boundaries by Superposition

Real problems give you all four edges. Decompose the prescribed data as

Problem T: top edge = $f_T(x)$ , other three = 0. We just solved this.
Problem B: bottom edge = $f_B(x)$ , other three = 0. Same formula with $y \to b - y$ .
Problem L: left edge = $f_L(y)$ , other three = 0. Same formula with $x \leftrightarrow y$ and $a \leftrightarrow b$ .
Problem R: right edge = $f_R(y)$ , other three = 0. Same as Problem L with the role of left/right swapped.

By linearity, the full solution is $u_T + u_B + u_L + u_R$ . Each piece is worked out independently using the single-boundary machinery we derived. There is no new mathematics — only careful labelling.

Why superposition is allowed

Laplace's equation is linear and homogeneous: $L[u_T + u_B + u_L + u_R] = L u_T + L u_B + L u_L + L u_R = 0$ . The boundary conditions add as well (each piece is non-zero on a different edge), so the sum exactly satisfies the full four-edge data. Uniqueness of the Dirichlet problem (proved in section 28-02) then guarantees this is the answer.

What the Decay Rates Really Mean

The factor $\sinh(n\pi y/a)/\sinh(n\pi b/a)$ is the secret carrier of every interesting physical fact about Laplace's equation in a rectangle. For large $n$ , both sinh's look like half their exponential, so

\dfrac{\sinh(n\pi y/a)}{\sinh(n\pi b/a)} \;\approx\; e^{-n\pi (b - y)/a}.

Read this carefully. At depth $(b - y)$ below the active top edge, the n-th boundary mode is suppressed by a factor of $e^{-n\pi (b - y)/a}$ . So a feature of characteristic width $\ell$ on the boundary (which has dominant frequency $n \sim a/\ell$ ) will decay over a depth of order $\ell/\pi$ . Narrow features die fast; broad features reach deep.

This is the same principle that lets you see a mountain from far away but not read a book on the mountaintop with binoculars: high spatial frequencies decay rapidly with distance.

Common Pitfalls

Wrong sign on the separation constant. If you accidentally set $X'' - \lambda X = 0$ (instead of $X'' + \lambda X = 0$ ) the $X$ solutions become hyperbolic, the side-wall conditions force $X \equiv 0$ , and the only solution you get is $u \equiv 0$ . Always let the geometry pick the sign for you.
Forgetting to divide by sinh(nπb/a). Without the normalisation, your formula equals $\sin(n\pi x/a) \cdot \sinh(n\pi b/a)$ on the top edge, which is huge for large $n$ . Your reconstruction of $f$ will look like a spike train.
Numerical overflow for tall rectangles. $\sinh(n\pi b/a)$ overflows around $n\pi b/a \approx 710$ in double precision. For a 4:1 aspect ratio that is around $n \approx 56$ . Either use the asymptotic $\sinh(z)/\sinh(z + \Delta) \approx e^{-\Delta}$ for large arguments, or truncate $N$ before overflow — the high-n terms are negligible anyway.
Boundary mismatch at the corners. If $f(0) \ne 0$ or $f(a) \ne 0$ , the basis cannot reproduce the corner values exactly and you get Gibbs ringing. Either subtract off a corner-linear interpolant before expanding, or live with the overshoot.
Using cosines instead of sines. Cosines are natural for Neumann conditions (zero normal derivative). For zero Dirichlet conditions on the side walls, only sines vanish at both endpoints — sines are forced. Picking the wrong basis ruins everything downstream.

Summary

The whole section in one breath:

Question	Answer
What PDE are we solving?	uₓₓ + u_yy = 0 on [0,a]×[0,b]
Which boundary is non-zero?	One edge at a time — superposition handles the rest
What ansatz?	Product u = X(x) Y(y)
What are the X-modes?	sin(nπx/a) for n = 1, 2, 3, ...
What are the Y-modes?	sinh(nπy/a) / sinh(nπb/a)
How are the coefficients chosen?	Fourier sine projection of f on the active edge
How fast do modes decay?	Like e^(-nπ(b−y)/a) — narrow features die fast
How do we implement it?	Compute b_n by quadrature, then one einsum

The single sentence to take with you

The Laplace problem on a rectangle is a Fourier-sine projection of the boundary along its homogeneous direction, paired with an exponential (sinh) interpolation across the inhomogeneous direction — geometry chooses the basis, orthogonality chooses the coefficients, and $\sinh$ chooses how deeply each boundary feature reaches into the interior.