Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

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

Explain why incompressible and irrotational flow forces both the velocity potential $\varphi$ and the stream function $\psi$ to satisfy Laplace's equation.
Read a flow picture: identify streamlines, equipotentials, stagnation points, and the dividing streamline that traces the body.
Use the complex potential $w(z) = \varphi + i\psi$ to combine elementary flows by simple addition.
Derive the flow past a circular cylinder from a uniform stream plus a doublet, and recover the surface pressure from Bernoulli.
Implement and verify a numerical Laplace solver and check it against the analytic answer.

Why this matters. Every clean equation in classical aerodynamics — lift on a thin wing, the Magnus force on a spinning ball, the shape of streamlines around a submarine hull — comes from one observation: an ideal fluid behaves like the gradient of a harmonic function. The same Laplace equation that solved electrostatics in the previous section solves fluid flow here, with the same techniques and the same intuitions.

The Big Picture

Imagine smoke pouring around a windshield, or water flowing past a bridge pier on a calm day. The fluid does something remarkable: it re-organizes itself into smooth, gently curving paths that hug the obstacle and then close back up behind it. There is no turbulence, no wake, no chaos — just a quiet, almost crystalline pattern of curves.

Look at any one of those curves and ask: what equation describes it? Astonishingly, the answer is the same equation that describes the steady-state temperature in a rod, the electrostatic potential between two charged plates, and the deflection of a soap film stretched across a wire loop:

\displaystyle \nabla^{2}\varphi \;=\; \frac{\partial^{2}\varphi}{\partial x^{2}} + \frac{\partial^{2}\varphi}{\partial y^{2}} \;=\; 0.

Different physics, different units, same equation. The reason is that Laplace's equation is the universal language of equilibrium in the absence of sources. Anywhere a quantity diffuses to its most balanced state, ∇² of that quantity is zero. Fluid flow is one more dialect of that language.

Intuition. A harmonic function is a function whose value at any point is the average of its values on any small circle around that point. Fluid in steady, source-free motion redistributes itself until it satisfies exactly that averaging property — locally, no point is special.

Two Assumptions: Incompressible & Irrotational

To get Laplace's equation we have to assume two things about the fluid. They are strong assumptions — real fluids violate both at high speeds — but they are exactly the right assumptions to expose the clean mathematics underneath.

1. Incompressible

The fluid's density is constant. No matter how the flow squeezes and stretches, every little parcel keeps the same volume. Mathematically this means the divergence of the velocity vector is zero:

\nabla \cdot \mathbf{v} \;=\; \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \;=\; 0.

Picture: draw a small square in the fluid. The amount flowing into the left side has to equal the amount flowing out of the right side (plus the same balance vertically), otherwise fluid would pile up. The divergence is just “net outflow per unit area” — and we are setting it to zero.

2. Irrotational

Tiny fluid parcels translate and stretch but do not spin. Insert a miniature pinwheel into the flow and it would not rotate. Mathematically the curl of the velocity is zero:

\nabla \times \mathbf{v} \;=\; \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} \;=\; 0.

Picture: the flow can curve and wind around obstacles all it wants — what matters is that locally there is no swirl. Fluid flowing around a cylinder still satisfies this: the curvature is global, but the local rotation is zero.

These two conditions are independent. Incompressibility is about volume; irrotationality is about local spin. A real river is roughly incompressible but very rotational (eddies everywhere). An idealised airflow over a slow-moving wing is approximately both.

The Velocity Potential φ

The curl-free condition $\nabla \times \mathbf{v} = 0$ has an immediate consequence from vector calculus: a curl-free field is the gradient of a scalar. Concretely, there exists a function $\varphi(x, y)$ — the velocity potential — such that

\mathbf{v} \;=\; \nabla \varphi, \qquad \text{i.e.}\quad u = \frac{\partial \varphi}{\partial x}, \;\; v = \frac{\partial \varphi}{\partial y}.

Why is this useful? Because instead of tracking two velocity components $u(x, y)$ and $v(x, y)$ we now only need to know one scalar function $\varphi$ . The flow has been compressed by a factor of two.

Analogy. Gravity is curl-free; its potential is height. Knowing the height map of a landscape tells you everything about where water will flow downhill. The velocity potential is the same idea, with the gradient pointing along the flow instead of perpendicular to contours of height.

The Stream Function ψ

Now use the other condition. Incompressibility $\nabla \cdot \mathbf{v} = 0$ says that $u$ and $v$ satisfy $\partial u / \partial x = -\partial v / \partial y$ . This is exactly the condition for the differential form $u\,dy - v\,dx$ to be exact, which means there is a function $\psi(x, y)$ — the stream function — with

u = \frac{\partial \psi}{\partial y}, \qquad v = -\frac{\partial \psi}{\partial x}.

The stream function has a magical geometric meaning. Compute the derivative of $\psi$ along a streamline (a curve that follows the velocity at every point):

\frac{d\psi}{dt} = \frac{\partial \psi}{\partial x}\frac{dx}{dt} + \frac{\partial \psi}{\partial y}\frac{dy}{dt} = -v\,u + u\,v = 0.

So $\psi$ is constant along streamlines. Streamlines are the level sets of ψ. If you want to draw the flow pattern, you just contour-plot $\psi$ .

Fluid never crosses a streamline. So any streamline can be replaced by a solid wall without changing the flow elsewhere. This is the trick that lets us solve flow past a cylinder using only point singularities: we find a superposition whose ψ = 0 streamline happens to be a circle, then declare that circle to be solid.

Why Both φ and ψ Satisfy Laplace's Equation

We have two scalar functions describing the flow. Each one inherits a differential equation from the other condition.

φ inherits incompressibility

Plug $u = \partial \varphi / \partial x$ and $v = \partial \varphi / \partial y$ into $\partial u / \partial x + \partial v / \partial y = 0$ :

\frac{\partial^{2}\varphi}{\partial x^{2}} + \frac{\partial^{2}\varphi}{\partial y^{2}} \;=\; \nabla^{2}\varphi \;=\; 0.

ψ inherits irrotationality

Plug $u = \partial \psi / \partial y$ and $v = -\partial \psi / \partial x$ into $\partial v / \partial x - \partial u / \partial y = 0$ :

-\frac{\partial^{2}\psi}{\partial x^{2}} - \frac{\partial^{2}\psi}{\partial y^{2}} \;=\; -\nabla^{2}\psi \;=\; 0.

So $\nabla^{2}\psi = 0$ too. Both functions are harmonic, and the entire arsenal from the earlier sections of this chapter — mean-value property, maximum principle, uniqueness with boundary conditions, separation of variables — applies directly.

Quantity	Defined by	PDE it satisfies	What its level sets are
Velocity potential φ	v = ∇φ (from curl = 0)	∇²φ = 0	Equipotential curves (⊥ to streamlines)
Stream function ψ	u = ∂ψ/∂y, v = -∂ψ/∂x (from div = 0)	∇²ψ = 0	Streamlines (along the flow)

Beautiful fact. The streamlines ψ = const and the equipotentials φ = const meet at right angles everywhere. This is a purely mathematical consequence of the Cauchy-Riemann equations below — and it is exactly what we see in the visualization above.

The Complex Potential w(z) (the Magic Trick)

Now we play a card that physicists love. Compare the two pairs of definitions:

u = \frac{\partial \varphi}{\partial x} = \frac{\partial \psi}{\partial y}, \qquad v = \frac{\partial \varphi}{\partial y} = -\frac{\partial \psi}{\partial x}.

Those are exactly the Cauchy-Riemann equations from complex analysis. The pair $(\varphi, \psi)$ behaves like the real and imaginary parts of a single complex-differentiable function of $z = x + iy$ . Define

w(z) \;=\; \varphi(x, y) + i\,\psi(x, y).

Then any analytic function of z is automatically a valid 2D potential flow. That is a stunning fact: complex analysis hands us an infinite catalog of pre-built solutions to Laplace's equation for free.

Even better: the velocity comes out by complex differentiation,

\frac{dw}{dz} \;=\; u - i\,v.

That sign on $v$ is what makes the magic work when you superimpose flows.

The catalogue of elementary complex potentials

Every interesting 2D flow in the rest of this section can be written as a sum of these atoms:

Uniform flow: $w(z) = U z$
Point source / sink: $w(z) = \frac{m}{2\pi}\log(z - z_0)$ (m > 0 = source, m < 0 = sink)
Point vortex: $w(z) = -i\,\frac{\Gamma}{2\pi}\log(z - z_0)$
Doublet: $w(z) = \frac{\mu}{2\pi}\cdot \frac{1}{z - z_0}$

Elementary Flows: The Building Blocks

Let's unpack each atom geometrically. For every one, ψ is constant along streamlines (so contour-plotting ψ gives you the flow pattern), and ∇²ψ = ∇²φ = 0 holds everywhere except at the singular point itself.

Flow	φ	ψ	Picture
Uniform (speed U, +x)	U x	U y	Parallel horizontal lines
Source (strength m)	(m/2π) ln r	(m/2π) θ	Rays radiating out
Sink (strength m)	-(m/2π) ln r	-(m/2π) θ	Rays radiating in
Vortex (circulation Γ)	(Γ/2π) θ	-(Γ/2π) ln r	Concentric circles
Doublet (strength μ)	(μ/2π) x/r²	-(μ/2π) y/r²	Cardioid-like loops

Notice the symmetry: for the source, the rays (streamlines) and the circles (equipotentials) are perpendicular — exactly as the Cauchy- Riemann equations demand. For the vortex, the roles are swapped: the circles are streamlines and the rays are equipotentials.

The doublet is the most subtle one. It is the limit of a source at

(-\epsilon, 0)

and a sink of equal strength at

(\epsilon, 0)

\epsilon \to 0

with the product

2 m \epsilon = \mu

held fixed. It is the 2D analogue of a dipole in electrostatics — and, as we will see, it is exactly what you add to a uniform flow to make a cylinder appear.

Superposition: Adding Flows Together

Because Laplace's equation is linear, the sum of two harmonic functions is harmonic. So the sum of two potential flows is itself a potential flow. Mechanically: if $w_1(z)$ and $w_2(z)$ are valid complex potentials, then so is $w_1(z) + w_2(z)$ .

This is the single most powerful idea in classical fluid mechanics. It means we can design flows by choosing what singularities to sprinkle into a uniform stream and watching what body shape emerges as the ψ = 0 streamline.

Think like an engineer. Want a flow that looks like wind passing a teardrop? Put a source upstream and a sink downstream inside a uniform flow; the closed streamlines form a Rankine oval. Want a circle? Replace the source/sink pair by a doublet. Want lift? Add a vortex. The geometry is encoded in the choice of singularities.

Interactive Flow Explorer

Try it. Move singularities around with the mouse, tune their strengths, and watch the streamlines (cyan) and equipotentials (yellow) reorganize. Notice that:

Streamlines never cross. A streamline can split at a stagnation point, but two distinct streamlines never share a non- stagnation point.
Streamlines and equipotentials are everywhere orthogonal — the Cauchy-Riemann equations in pixels.
Stagnation points (where v = 0) appear as small yellow crosshairs. They are the points where the dividing streamline meets the body.
In Flow past cylinder, the ψ = 0 streamline coincides exactly with the cylinder surface — that is why the cylinder is a valid “wall”.
In Lifting cylinder (Magnus), adding a vortex pushes the stagnation points around the cylinder. Asymmetry between top and bottom is exactly what creates lift in a spinning ball — the same physics behind a curving free kick.

Potential-Flow Explorer

Drag any marker. Cyan = streamlines (ψ = const). Yellow = equipotentials (φ = const).

Preset

Uniform flow U1.00

Cylinder radius a1.00

Contour density28

Add

doublet(0.00, 0.00)

strength6.28

Tip: drag a marker to move it. Each elementary flow individually satisfies Laplace's equation; their sum still satisfies Laplace's equation because $\nabla^{2}$ is linear. The cyan curves are level sets of $\psi$ ; fluid flows along them and never across.

Worked Example: Flow Past a Cylinder

Let's do the most famous superposition in fluid mechanics from scratch, by hand. The goal: derive the flow around a circular cylinder of radius $a$ sitting in a uniform stream of speed $U$ , and find where on the surface the pressure is highest and lowest.

▸ Click to expand the full hand-derivation

Step 1: Pick the singularities. A uniform flow alone gives flat horizontal streamlines — no body. We need to distort those lines into a closed circle. The cheapest singularity with the right symmetry (no net source, no preferred direction other than along x) is a doublet at the origin.

Step 2: Write the candidate stream function.

\psi(x, y) \;=\; U y \;+\; \Big(-\frac{\mu}{2\pi}\Big)\frac{y}{x^{2} + y^{2}} \;=\; U y\Big(1 - \frac{\mu / (2 \pi U)}{x^{2} + y^{2}}\Big).

Define $a^{2} = \mu / (2 \pi U)$ , equivalently $\mu = 2 \pi U a^{2}$ . Then

\psi \;=\; U y \Big(1 - \frac{a^{2}}{r^{2}}\Big), \quad r^{2} = x^{2} + y^{2}.

Step 3: Find the ψ = 0 streamline. Setting $\psi = 0$ gives $y = 0$ (the x-axis) or $r^{2} = a^{2}$ (the circle of radius a). That circle is a valid solid wall — we just declare it to be the cylinder.

Step 4: Velocity on the surface. In polar coordinates the streamline expression becomes $\psi = U r \sin\theta (1 - a^{2}/r^{2})$ . Compute the tangential velocity on the cylinder $r = a$ :

v_{\theta} = -\frac{\partial \psi}{\partial r}\bigg|_{r=a} = -U \sin\theta\Big(1 + \frac{a^{2}}{r^{2}}\Big)\bigg|_{r=a} = -2U \sin\theta.

Step 5: Stagnation points. $v_{\theta} = 0$ at $\theta = 0$ and $\theta = \pi$ — the front and rear of the cylinder. Symmetric, as you would expect from a symmetric flow.

Step 6: Pressure from Bernoulli. For steady incompressible flow, $p + \tfrac{1}{2}\rho |\mathbf{v}|^{2} = \text{const}$ . With the far-field reference $p_\infty, U$ :

p(\theta) - p_{\infty} = \tfrac{1}{2} \rho \big(U^{2} - v_{\theta}^{2}\big) = \tfrac{1}{2} \rho U^{2} \big(1 - 4 \sin^{2}\theta\big).

So pressure is highest at the front and rear stagnation points (θ = 0 and π), and lowest at the sides (θ = π/2 and 3π/2), where the surface speed reaches $2U$ .

Step 7: Net force on the cylinder. Integrate $p(\theta)$ times the surface normal around the cylinder. Both the drag (x-component) and the lift (y-component) integrate to zero by symmetry. Net force = 0. This is d'Alembert's paradox.

Sanity check the strength formula numerically: $\mu = 2\pi \cdot 1 \cdot 1^{2} \approx 6.2832$ — exactly the doublet strength used in the “Flow past cylinder” preset above.

Plain Python: Building Streamlines from Scratch

Before any PyTorch, let's implement the cylinder flow with nothing but NumPy and a contour plot. Reading this code line-by-line cements the math: each formula in the worked example is a single line of vectorized array arithmetic.

Flow past a cylinder via uniform flow + doublet

🐍cylinder_flow.py

Explanation(23)

Code(45)

1Import NumPy

NumPy gives us vectorized array math. Every quantity below (ψ, φ, velocity) is computed on the whole grid in one shot, no Python loop.

2Import matplotlib

We need contour plots: every streamline is a level set ψ = const, and every equipotential is a level set φ = const. `contour` draws those level curves directly.

5Grid resolution

400 × 250 grid points in world coordinates. Higher resolution gives smoother contours but more memory; this is plenty for a single illustrative figure.

EXECUTION STATE

nx = 400

ny = 250

6x-axis samples

We span x ∈ [-4, 4]. The cylinder sits at the origin, so this gives ~3 cylinder radii of free flow on either side — enough to see how streamlines reform downstream.

EXAMPLE

x[0] = -4.0, x[-1] = 4.0

7y-axis samples

y ∈ [-2.5, 2.5]. We keep the aspect ratio close to the canvas so the cylinder looks round.

8Build the 2D coordinate arrays

`meshgrid` turns the 1-D x and y vectors into two 2-D arrays X[j,i], Y[j,i] holding the (x, y) coordinates of every grid point. After this line we can write any function of (x, y) once and get a whole 2-D array back.

EXAMPLE

X.shape == (250, 400),  Y.shape == (250, 400)

11Free-stream speed

The fluid at infinity moves with horizontal speed U. Everything downstream is measured relative to this.

12Cylinder radius

We will discover that the choice μ = 2π U a² in line 13 makes the dividing streamline ψ = 0 coincide with the circle r = a. Pick a, then μ falls out.

13Doublet strength

A doublet is a source and a sink pushed infinitesimally close together. Its stream function ψ = -(μ / 2π) · y/r² has a pole at the origin. The magic: when you add it to U·y you get ψ_total = U y (1 - a²/r²), and ψ_total = 0 on the circle r = a.

EXAMPLE

μ = 2π · 1.0 · 1.0² ≈ 6.2832

16Squared radius

We need r² in many places (denominators for doublet, source, vortex). Computing it once and reusing is faster than calling np.hypot repeatedly.

17Mask the interior

Points inside the cylinder are not part of the physical fluid domain. We mark them so the contour routine knows to stop the lines at the surface.

20Uniform-flow stream function

For flow at speed U in the +x direction, ψ = U y. Lines of constant y are horizontal streamlines — exactly what you expect for uniform flow.

EXAMPLE

ψ_uniform = [[-2.5, ..., -2.5], ..., [2.5, ..., 2.5]] (constant along each row)

21Doublet stream function

ψ_doublet(x, y) = -(μ / 2π) · y / r². This expression vanishes on the x-axis and goes to infinity at the origin. It is the only solution of Laplace's equation that decays like 1/r and has dipole symmetry.

22Superposition

Add them. Because Laplace's equation is linear (∇²(ψ₁ + ψ₂) = ∇²ψ₁ + ∇²ψ₂), the sum still solves Laplace's equation. This is the single most powerful idea in potential flow: complex flows are sums of elementary ones.

23Erase the interior

Setting masked entries to NaN tells matplotlib's `contour` to skip them, so streamlines stop exactly at the cylinder.

26Uniform-flow potential

Velocity potential satisfies u = ∂φ/∂x, v = ∂φ/∂y. For uniform horizontal flow u = U, v = 0, so φ = U x. Equipotentials are vertical lines x = const.

27Doublet potential

φ_doublet(x, y) = (μ / 2π) · x / r². Notice the swap: doublet has y/r² for ψ and x/r² for φ. This swap is the Cauchy-Riemann equations in disguise — see the worked example.

32Verify Laplace's equation

We approximate ∇²ψ by the standard 5-point stencil and check that, away from the origin, it is essentially zero. If our formulas were wrong, this would scream at us with values of order 1 instead of 1e-12.

EXAMPLE

max |Laplacian psi| ~ 1e-3 in finite-precision

36Print check

We slice off the first/last 10 rows and columns so the periodic-wrap from np.roll doesn't contaminate the result, and we skip a band near the origin where the doublet pole produces large numerical error.

41Plot streamlines

`contour` draws 25 level sets of ψ. These are the streamlines — fluid particles travel along them and never cross them (because ψ is conserved along the flow).

42Plot equipotentials

25 level sets of φ. Where streamlines and equipotentials cross, they cross at 90° — that is the orthogonality property of conjugate harmonic functions.

44Draw the cylinder body

A pink translucent disk so the eye immediately reads the geometry. The cylinder boundary is itself the streamline ψ = 0.

45Equal aspect ratio

Without this, a perfect circle would render as an ellipse and the perpendicularity of ψ and φ contours would be lost.

22 lines without explanation

1import numpy as np
2import matplotlib.pyplot as plt
3
4# --- Build a grid in world coordinates ---
5nx, ny = 400, 250
6x = np.linspace(-4, 4, nx)
7y = np.linspace(-2.5, 2.5, ny)
8X, Y = np.meshgrid(x, y)
9
10# --- Elementary flow: uniform stream + doublet (cylinder of radius a) ---
11U  = 1.0       # free-stream speed
12a  = 1.0       # cylinder radius
13mu = 2 * np.pi * U * a ** 2  # doublet strength so the cylinder appears
14
15# Mask the inside of the cylinder so contours stop at the surface.
16R2   = X ** 2 + Y ** 2
17mask = R2 < a ** 2
18
19# --- Stream function psi(x, y) ---
20psi_uniform = U * Y
21psi_doublet = -mu / (2 * np.pi) * Y / R2
22psi         = psi_uniform + psi_doublet
23psi[mask]   = np.nan
24
25# --- Velocity potential phi(x, y) ---
26phi_uniform = U * X
27phi_doublet = mu / (2 * np.pi) * X / R2
28phi         = phi_uniform + phi_doublet
29phi[mask]   = np.nan
30
31# --- Check Laplace: nabla^2 psi should be ~0 away from origin ---
32lap = (np.roll(psi, -1, 0) + np.roll(psi, 1, 0)
33       + np.roll(psi, -1, 1) + np.roll(psi, 1, 1)
34       - 4 * psi) / ((x[1] - x[0]) ** 2)
35print("max |Laplacian psi| (interior, away from origin):",
36      np.nanmax(np.abs(lap[10:-10, 10:-10])))
37
38# --- Plot streamlines and equipotentials ---
39fig, ax = plt.subplots(figsize=(8, 5))
40ax.contour(X, Y, psi, levels=25, colors="#38bdf8", linewidths=0.8)
41ax.contour(X, Y, phi, levels=25, colors="#fbbf24",
42           linewidths=0.5, alpha=0.6)
43ax.add_patch(plt.Circle((0, 0), a, color="#f472b6", alpha=0.25))
44ax.set_aspect("equal")
45plt.show()

What you should see when you run it

The printed Laplacian residual is of order $10^{-3}$ (limited by the 5-point stencil discretization, not by our formulas). The plot shows cyan streamlines flowing around a pink cylinder, yellow equipotentials crossing them at right angles, and front/rear stagnation points where the dividing streamline meets the surface.

PyTorch: Solving Laplace by Relaxation

The closed-form is beautiful — but in practice we will not always have one. The general recipe is to discretize $\nabla^{2} \varphi = 0$ on a grid, fix Dirichlet boundary data, and iterate. Jacobi relaxation is the textbook starting point: replace every interior value with the average of its four neighbours and repeat until nothing changes.

PyTorch makes this trivial because the four neighbour shifts are $\texttt{torch.roll}$ calls, and the averaging is a single fused arithmetic op. The exact same code runs on a GPU.

Jacobi relaxation for Laplace's equation around a cylinder

🐍laplace_relaxation.py

Explanation(21)

Code(62)

1Import torch

We use PyTorch tensors (not autograd) here. They give us GPU support and vectorized neighbour shifts via `torch.roll`. The same code runs on CPU and CUDA by flipping `device`.

3Function signature

Solve Laplace's equation on a (nx, ny) grid using Jacobi relaxation. We compare the result against the analytic flow-past-cylinder potential as a numerical sanity check.

16x and y axes

Same box as the explorer above. `linspace` gives evenly-spaced sample points including the endpoints, so boundary nodes are exactly on x = ±4 and y = ±2.5.

18Build the meshgrid

`indexing="xy"` gives X.shape == (ny, nx), matching standard Cartesian convention. Now X[j, i] is the x-coordinate and Y[j, i] is the y-coordinate of the grid point (j, i).

EXECUTION STATE

X.shape = (120, 200)

Y.shape = (120, 200)

22Cylinder radius

We solve flow past a unit cylinder. Picking a = 1 keeps the numbers clean.

23Build the interior mask

Boolean tensor that is True inside the cylinder. We will (a) force φ = 0 there and (b) refuse to update those nodes in the iterations.

26Analytic potential

Closed-form answer derived in the worked example: φ = U x + (μ / 2π) · x / r². We use it both as the Dirichlet boundary data and as the ground-truth to measure error.

EXAMPLE

At (x, y) = (2, 0):  φ = 2 + 0.5 = 2.5

32Initial guess

Start φ at zero everywhere. Jacobi relaxation will spread the boundary values inward until they satisfy ∇²φ = 0. The initial guess is irrelevant — Laplace's equation has a unique solution given Dirichlet boundary data.

34Top / bottom edges

Copy the analytic potential on rows 0 and -1. These rows will not be touched by the iteration.

36Left / right edges

Copy the analytic potential on the leftmost and rightmost columns. With all four edges fixed, the boundary-value problem is well-posed and has a unique interior solution.

38Zero inside the cylinder

We treat the interior of the cylinder as a hole. The cylinder surface itself is part of the boundary (carries the analytic value, which equals 0 on the circle for this superposition).

41Build the interior mask

`interior` is True at every node we are allowed to update. We start with all-True and switch off the 4 edges plus the cylinder mask.

47Jacobi update — definition

Jacobi's rule is brutally simple: at each step replace every interior φ value with the average of its 4 neighbours. That average is exactly the 5-point discrete Laplacian set to zero, so a fixed point of the iteration solves ∇²φ = 0 to discretization order.

48Shift north

`torch.roll(phi, -1, 0)` shifts every row up by one, so position (j, i) now holds the value that used to live at (j+1, i) — the north neighbour. Edge wraparound is harmless because we will overwrite the boundary rows.

49Shift south

`+1` on axis 0 shifts down: every cell sees its south neighbour.

50Shift east

Axis 1 is the column axis. `-1` shifts left, so position (j, i) sees its east neighbour.

51Shift west

`+1` on axis 1: every cell sees its west neighbour.

52Average the 4 neighbours

Vectorized 5-point Laplace solve in a single tensor op. No Python loop over (i, j). On a GPU this is essentially free.

EXAMPLE

phi_new[j, i] = 0.25 * (phi[j-1, i] + phi[j+1, i] + phi[j, i+1] + phi[j, i-1])

53Apply only on the interior

`torch.where(interior, phi_new, phi)` keeps the new value where we are allowed to update, and preserves the old (boundary / cylinder) value elsewhere. This is how we honour the Dirichlet boundary.

55Periodic error report

Every 2000 iterations we print the mean absolute deviation from the analytic potential. After ~8000 iterations on a 200 × 120 grid it drops below 1e-3, which is the floor set by our discretization (∼ Δx²).

56Compute the error tensor

Element-wise difference, absolute value, restricted to the interior, then mean. We do not compare on the boundary because we hardcoded the boundary to be exactly correct.

41 lines without explanation

1import torch
2
3def solve_laplace_2d(
4    nx: int = 200,
5    ny: int = 120,
6    iters: int = 8000,
7    device: str = "cpu",
8) -> torch.Tensor:
9    """
10    Solve nabla^2 phi = 0 on a rectangular box with Dirichlet
11    boundary conditions chosen so the analytic answer is uniform
12    flow + doublet (i.e. flow past a cylinder of radius a = 1).
13    We then check phi against the analytic potential.
14    """
15    # World box: x in [-4, 4], y in [-2.5, 2.5]
16    x = torch.linspace(-4, 4, nx, device=device)
17    y = torch.linspace(-2.5, 2.5, ny, device=device)
18    X, Y = torch.meshgrid(x, y, indexing="xy")
19    R2   = X ** 2 + Y ** 2
20
21    # Cylinder mask (interior is excluded from the solution).
22    a    = 1.0
23    mask = R2 < a ** 2
24
25    # Analytic potential phi = U*x + (mu / 2pi) * x / r^2  (with U = 1).
26    U   = 1.0
27    mu  = 2 * torch.pi * U * a ** 2
28    phi_exact = U * X + (mu / (2 * torch.pi)) * X / R2
29    phi_exact[mask] = 0.0
30
31    # --- Set up the iterative (Jacobi) solver ---
32    phi = torch.zeros_like(X)
33    # Boundary condition: copy the analytic answer on the outer edge
34    # and on the cylinder surface. The interior must be discovered.
35    phi[0, :]  = phi_exact[0, :]
36    phi[-1, :] = phi_exact[-1, :]
37    phi[:, 0]  = phi_exact[:, 0]
38    phi[:, -1] = phi_exact[:, -1]
39    phi[mask]  = 0.0  # treat cylinder interior as a "hole"
40
41    # Build an interior mask (everything we are allowed to update).
42    interior = torch.ones_like(X, dtype=torch.bool)
43    interior[0, :] = interior[-1, :] = False
44    interior[:, 0] = interior[:, -1] = False
45    interior[mask] = False
46
47    # --- Jacobi iteration: phi_new = average of 4 neighbours ---
48    for k in range(iters):
49        north = torch.roll(phi, -1, 0)
50        south = torch.roll(phi,  1, 0)
51        east  = torch.roll(phi, -1, 1)
52        west  = torch.roll(phi,  1, 1)
53        phi_new = 0.25 * (north + south + east + west)
54        phi = torch.where(interior, phi_new, phi)
55
56        if (k + 1) % 2000 == 0:
57            err = (phi - phi_exact)[interior].abs().mean()
58            print(f"iter {k+1:5d}   mean |phi - phi_exact| = {err.item():.4e}")
59
60    return phi
61
62phi = solve_laplace_2d()

What you should see in your terminal

After 8000 Jacobi sweeps on a 200 × 120 grid the mean error settles near $5 \times 10^{-3}$ — a finite-difference discretization-floor that goes away with finer grids or smarter solvers (successive over-relaxation, conjugate gradient, multigrid). The numerical φ matches the analytic flow-past-cylinder potential everywhere except in a thin numerical boundary layer near the cylinder surface.

d'Alembert's Paradox

Integrate the surface pressure around the cylinder. By symmetry the net force is zero: no drag, no lift. Yet anyone who has stuck a hand out a car window knows that flow past a cylinder absolutely produces drag.

This contradiction, first noticed by d'Alembert in 1752, is not a bug in our derivation. It is a precise statement of where ideal potential flow stops being a faithful model of reality:

Real fluids have viscosity. Even a tiny amount creates a thin boundary layer at the surface where rotation builds up. Vorticity gets shed downstream as a wake.
The wake breaks the symmetry. Pressure on the rear of the cylinder drops below pressure on the front, and the integral no longer cancels.
So why study potential flow? Because outside the boundary layer and the wake, real flows look exactly like the potential prediction. Aircraft wings, propellers, hydrofoils and many sonar designs all start their life in the language we just built.

The next chapter on the Navier-Stokes equations is where viscosity comes back in and the paradox dissolves.

Where This Shows Up

Field	How Laplace + potential flow is used
Aerodynamics	Thin-airfoil theory and the Kutta-Joukowski lift theorem use the complex potential to predict lift on wings before any CFD code runs.
Naval architecture	Submarine and ship-hull design uses panel methods — distributing thousands of source/sink panels on the hull and solving for their strengths to match the surface boundary condition.
Groundwater hydrology	Darcy's law makes pressure in a porous aquifer satisfy ∇²p = 0. Source = well drawing water out; sink = recharge basin. Same math, decades-long timescales.
Sports physics	The Magnus force on a spinning football or tennis ball is exactly the lifting-cylinder solution in 2D — visible above when you toggle the vortex strength on top of the cylinder doublet.
Computer graphics	Smooth, divergence-free vector fields for fluid simulation in animation are often initialized from harmonic potentials and then projected back onto the divergence-free subspace.
Machine learning	Physics-informed neural networks (PINNs) learn φ(x, y) by penalizing the Laplacian residual at sampled points. The very same loss function generalizes from electrostatics to fluid flow with no change of code.

The takeaway. Once you can recognize that some physical quantity is the gradient of a scalar that diffuses to its most-balanced state, you do not need to re-derive anything. Reach for Laplace's equation, pick a boundary condition, and apply whichever technique from this chapter fits the geometry — separation of variables in a box, in a circle, or numerical relaxation everywhere else.

Summary

For an incompressible & irrotational 2D flow, both the velocity potential $\varphi$ and the stream function $\psi$ are harmonic: $\nabla^{2}\varphi = \nabla^{2}\psi = 0$ .
Streamlines = level sets of $\psi$ . Equipotentials = level sets of $\varphi$ . They intersect at right angles (Cauchy-Riemann).
The complex potential $w(z) = \varphi + i\psi$ is analytic. Every analytic function is a valid 2D potential flow, so complex analysis gives us an infinite library of solutions.
Elementary atoms — uniform flow, source/sink, vortex, doublet — can be superposed (added) because Laplace's equation is linear. Body shapes appear automatically as the $\psi = 0$ streamlines.
Uniform flow + doublet at the origin gives the famous flow past a cylinder. Adding a vortex tilts the stagnation points and produces the Magnus lift.
Numerically, Laplace's equation surrenders to Jacobi relaxation: replace every interior value with the average of its four neighbours and iterate. PyTorch makes this a three-line loop that runs on GPU.
The theory predicts zero drag (d'Alembert's paradox); that is the precise place where viscosity must enter, leading us into the Navier-Stokes equations of the next chapter.