Module 2.1
Partial differential equations are the mathematical language of continuum physics. Every model of heat conduction, sound propagation, fluid flow, electromagnetic radiation, and quantum mechanics ultimately reduces to a PDE. Classifying a PDE as elliptic, parabolic, or hyperbolic is not merely a taxonomic exercise: the classification determines which boundary and initial conditions produce a well-posed problem, which solution techniques apply, and which numerical schemes will converge. Mastering this classification and the canonical solution strategies (separation of variables, Fourier series, d'Alembert's formula) gives you a toolkit that transfers directly to every applied science and engineering discipline. Furthermore, understanding well-posedness in the sense of Hadamard guards against wasting effort on ill-posed formulations that no algorithm can salvage.
After completing this module you will be able to:
A second-order linear PDE in two independent variables \(x,y\) is an equation of the form
$$ A\,u_{xx} + 2B\,u_{xy} + C\,u_{yy} + D\,u_x + E\,u_y + F\,u = G, $$where \(A,B,C,D,E,F,G\) may depend on \(x\) and \(y\) but not on \(u\) or its derivatives.
The discriminant of the principal part is \(\Delta = B^2 - AC\). The PDE is classified at each point as:
A twice continuously differentiable function \(u : \Omega \to \mathbb{R}\) is harmonic on the open set \(\Omega\) if it satisfies Laplace's equation everywhere in \(\Omega\):
$$ \nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0. $$A problem consisting of a PDE together with boundary and/or initial conditions is well-posed in the sense of Hadamard if:
Let \(f : [0,L] \to \mathbb{R}\) be piecewise smooth. Its Fourier sine series is
$$ f(x) \sim \sum_{n=1}^{\infty} b_n \sin\!\left(\frac{n\pi x}{L}\right), \qquad b_n = \frac{2}{L}\int_0^L f(x)\sin\!\left(\frac{n\pi x}{L}\right)dx. $$The separation of variables technique assumes a product solution \(u(x,t) = X(x)\,T(t)\), substitutes into the PDE, and separates the resulting expression into two ODEs, each equal to a constant \(-\lambda\) called the separation constant.
Let \(\Omega \subset \mathbb{R}^n\) be a bounded, connected open set and let \(u \in C^2(\Omega) \cap C(\overline{\Omega})\) be harmonic in \(\Omega\). Then
$$ \max_{\overline{\Omega}} u = \max_{\partial\Omega} u. $$In particular, a non-constant harmonic function cannot attain its maximum in the interior of \(\Omega\).
Proof sketch. Suppose \(u\) attains its maximum at an interior point \(x_0 \in \Omega\). By the mean-value property of harmonic functions, \(u(x_0)\) equals the average of \(u\) over any ball \(B_r(x_0) \subset \Omega\). Since \(u(x_0)\) is the maximum, every value in the ball must equal \(u(x_0)\). By connectedness, \(u\) is constant on \(\Omega\), contradicting the hypothesis.
Consider the heat equation \(u_t = k\,u_{xx}\) on \((0,L)\times(0,\infty)\) with initial condition \(u(x,0)=f(x)\) and homogeneous Dirichlet boundary conditions \(u(0,t)=u(L,t)=0\). If \(u_1\) and \(u_2\) are two classical solutions, then \(u_1 \equiv u_2\).
Proof. Let \(w = u_1 - u_2\). Then \(w\) satisfies \(w_t = k\,w_{xx}\) with \(w(x,0)=0\) and \(w(0,t)=w(L,t)=0\). Define the energy
$$ E(t) = \frac{1}{2}\int_0^L w(x,t)^2\,dx. $$Then \(E(0)=0\) and
$$ E'(t) = \int_0^L w\,w_t\,dx = k\int_0^L w\,w_{xx}\,dx = -k\int_0^L w_x^2\,dx \leq 0, $$where we integrated by parts and used the boundary conditions. Since \(E(t)\geq 0\) and \(E'(t)\leq 0\) with \(E(0)=0\), we conclude \(E(t)=0\) for all \(t\geq 0\), hence \(w\equiv 0\).
For equations with variable coefficients, the type can change from point to point. The Tricomi equation \(y\,u_{xx}+u_{yy}=0\) is elliptic for \(y>0\), parabolic on \(y=0\), and hyperbolic for \(y<0\).
Separation of variables works only for special geometries and boundary conditions. Most PDEs arising in practice must be solved numerically. The technique in this module provides insight and benchmarks, not a universal solver.
At discontinuities, the Fourier series converges to the midpoint of the jump and exhibits the Gibbs phenomenon: oscillations of roughly 9% overshoot that do not disappear as more terms are added, although they become narrower.
The classical wave equation \(u_{tt}=c^2 u_{xx}\) assumes small-amplitude vibrations of a perfectly flexible, uniform string. Real vibrating beams obey the biharmonic equation \(u_{tt}+\gamma^2 u_{xxxx}=0\); real membranes involve the 2D wave equation with tension dependent on geometry.
Well-posedness (Hadamard) is a necessary condition for a good numerical scheme, but not sufficient. A well-posed problem may still be ill-conditioned, so that floating-point perturbations are amplified. Stability analysis of the numerical scheme is a separate requirement (see Module 2.2).
Problem. Solve the initial-boundary value problem
$$ u_t = k\,u_{xx}, \quad 0 < x < L,\; t > 0, $$ $$ u(0,t)=0,\quad u(L,t)=0,\quad u(x,0)=f(x)=x(L-x). $$Step 1 — Separation of variables. Assume \(u(x,t)=X(x)T(t)\). Substituting:
$$ X(x)T'(t) = k\,X''(x)T(t) \implies \frac{T'}{kT} = \frac{X''}{X} = -\lambda. $$Step 2 — Spatial ODE. \(X'' + \lambda X = 0\) with \(X(0)=X(L)=0\). Non-trivial solutions require \(\lambda_n = (n\pi/L)^2\) with eigenfunctions \(X_n(x) = \sin(n\pi x/L)\), \(n=1,2,3,\ldots\)
Step 3 — Temporal ODE. \(T' + k\lambda_n T = 0 \implies T_n(t) = e^{-k(n\pi/L)^2 t}\).
Step 4 — Superposition.
$$ u(x,t) = \sum_{n=1}^{\infty} b_n \sin\!\left(\frac{n\pi x}{L}\right) e^{-k(n\pi/L)^2 t}. $$Step 5 — Fourier coefficients.
$$ b_n = \frac{2}{L}\int_0^L x(L-x)\sin\!\left(\frac{n\pi x}{L}\right)dx. $$Integration by parts (twice) gives:
$$ b_n = \frac{4L^2}{n^3\pi^3}\bigl[1-(-1)^n\bigr] = \begin{cases} \dfrac{8L^2}{n^3\pi^3} & n \text{ odd},\\[4pt] 0 & n \text{ even}.\end{cases} $$Final solution:
$$ u(x,t) = \frac{8L^2}{\pi^3}\sum_{\substack{n=1,3,5,\ldots}}^{\infty} \frac{1}{n^3}\sin\!\left(\frac{n\pi x}{L}\right)e^{-k(n\pi/L)^2 t}. $$Problem. Solve
$$ u_{tt} = c^2\,u_{xx}, \quad -\infty < x < \infty,\; t > 0, $$ $$ u(x,0) = \phi(x), \quad u_t(x,0)=\psi(x). $$Step 1 — Change of variables. Let \(\xi = x-ct\) and \(\eta = x+ct\). The PDE becomes \(u_{\xi\eta}=0\), whose general solution is \(u = F(\xi)+G(\eta) = F(x-ct)+G(x+ct)\).
Step 2 — Apply initial conditions.
$$ F(x)+G(x) = \phi(x), $$ $$ -c\,F'(x)+c\,G'(x) = \psi(x). $$Step 3 — Integrate the second equation:
$$ G(x)-F(x) = \frac{1}{c}\int_0^x \psi(s)\,ds + K. $$Step 4 — Solve for \(F\) and \(G\):
$$ F(x) = \frac{1}{2}\phi(x) - \frac{1}{2c}\int_0^x \psi(s)\,ds - \frac{K}{2}, $$ $$ G(x) = \frac{1}{2}\phi(x) + \frac{1}{2c}\int_0^x \psi(s)\,ds + \frac{K}{2}. $$D'Alembert's formula:
$$ u(x,t) = \frac{1}{2}\bigl[\phi(x-ct)+\phi(x+ct)\bigr] + \frac{1}{2c}\int_{x-ct}^{x+ct}\psi(s)\,ds. $$The first term represents two copies of the initial displacement traveling in opposite directions at speed \(c\). The second term accumulates the effect of the initial velocity over the domain of dependence \([x-ct,\,x+ct]\).
Problem. Classify the PDE \(3u_{xx}+4u_{xy}+u_{yy}-2u_x+u=0\).
Step 1 — Identify coefficients: \(A=3\), \(2B=4\) so \(B=2\), \(C=1\).
Step 2 — Compute discriminant: \(\Delta=B^2-AC=4-3=1>0\).
Conclusion: The equation is hyperbolic at every point (since the coefficients are constant).
The square wave on \([-\pi,\pi]\) defined by \(f(x)=\begin{cases}1&0<x<\pi\\-1&-\pi<x<0\end{cases}\) has the Fourier series
$$ f(x) \sim \frac{4}{\pi}\sum_{k=0}^{\infty}\frac{1}{2k+1}\sin\bigl((2k+1)x\bigr). $$The code below plots partial sums for increasing \(N\) so you can observe convergence and the Gibbs phenomenon near the discontinuities.
Output will display four plots showing Fourier partial sums with 1, 3, 9, and 49 odd terms converging to the square wave.
We frame PDE classification as a reinforcement-learning problem. An autonomous agent observes physical phenomena and must select the correct mathematical model and boundary conditions.
The agent's state is a feature vector describing the observed physical system:
The agent selects from a discrete action space:
Over many episodes the agent discovers the mapping between physical features and PDE type. It learns, for instance, that the presence of a single first-order time derivative with a second-order spatial derivative signals a parabolic equation, while two time derivatives signal a hyperbolic one. It also learns that elliptic equations require boundary data on the entire boundary (no initial conditions), while parabolic equations need an initial condition plus boundary data, and hyperbolic equations need initial data for both \(u\) and \(u_t\). The Q-table converges to a policy that simultaneously classifies, poses, and solves PDEs with high reliability.
Exercise 2.1.1. Classify each PDE as elliptic, parabolic, or hyperbolic:
(a) \(A=1,B=0,C=1\). \(\Delta=0-1=-1<0\). Elliptic.
(b) \(A=1,B=0,C=-1\). \(\Delta=0+1=1>0\). Hyperbolic.
(c) \(A=1,B=1,C=1\). \(\Delta=1-1=0\). Parabolic.
(d) At \((1,2)\): \(A=4,B=-2,C=1\). \(\Delta=4-4=0\). Parabolic. (In fact, \(\Delta=0\) for all \((x,y)\neq(0,0)\).)
Exercise 2.1.2. Solve the heat equation \(u_t = u_{xx}\) on \(0
By linearity and the eigenfunctions of \(-\partial_{xx}\) with Dirichlet BCs:
$$ u(x,t) = e^{-t}\sin(x) + 3e^{-4t}\sin(2x). $$The first mode decays with rate 1, the second with rate 4. Higher-frequency modes decay faster.
Exercise 2.1.3. Use d'Alembert's formula to solve \(u_{tt}=4u_{xx}\) on \(\mathbb{R}\) with \(u(x,0)=e^{-x^2}\) and \(u_t(x,0)=0\).
Here \(c=2\), \(\phi(x)=e^{-x^2}\), \(\psi(x)=0\). D'Alembert gives
$$ u(x,t) = \frac{1}{2}\left[e^{-(x-2t)^2}+e^{-(x+2t)^2}\right]. $$Two Gaussian pulses move in opposite directions at speed 2.
Exercise 2.1.4. Prove that if \(u\) is harmonic in the disk \(x^2+y^2<1\) and continuous on the closed disk with \(u=\cos^2\theta\) on the boundary, then \(\frac{1}{2}\leq u(x,y)\leq 1\) inside the disk. (Hint: use \(\cos^2\theta=\frac{1}{2}+\frac{1}{2}\cos 2\theta\) and the maximum principle.)
On the boundary, \(\min \cos^2\theta = 0\) and \(\max \cos^2\theta = 1\). By the maximum principle, \(0 \leq u \leq 1\) in the closed disk. For a sharper bound: the mean-value property gives \(u(0,0) = \frac{1}{2\pi}\int_0^{2\pi}\cos^2\theta\,d\theta = \frac{1}{2}\). Since \(u\) is harmonic and non-constant, the strong minimum principle ensures \(u > 0\) in the interior. The explicit Poisson-integral solution gives \(u(r,\theta)=\frac{1}{2}+\frac{r^2}{2}\cos 2\theta\), confirming \(\frac{1}{2}-\frac{r^2}{2}\leq u \leq \frac{1}{2}+\frac{r^2}{2}\leq 1\).
Exercise 2.1.5. Write a Python function that computes the Fourier sine coefficients \(b_n\) numerically (via scipy.integrate.quad) for an arbitrary function \(f\) on \([0,L]\) and returns the partial sum \(S_N(x)\).
Exercise 2.1.6. Implement the analytic solution of the heat equation from Example 2.1.1 and produce an animation showing how \(u(x,t)\) decays in time for \(L=\pi\), \(k=0.5\). Use at least 20 terms in the series.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
L, k_coeff = np.pi, 0.5
x = np.linspace(0, L, 300)
N_terms = 50
def u_heat(x, t):
s = np.zeros_like(x)
for m in range(N_terms):
n = 2 * m + 1
bn = 8 * L**2 / (n**3 * np.pi**3)
s += bn * np.sin(n * np.pi * x / L) * np.exp(-k_coeff * (n * np.pi / L)**2 * t)
return s
fig, ax = plt.subplots()
line, = ax.plot(x, u_heat(x, 0), 'b-')
ax.set_ylim(-0.1, L**2/4 + 0.2)
ax.set_xlabel('x'); ax.set_ylabel('u(x,t)')
title = ax.set_title('t = 0.00')
def animate(frame):
t = frame * 0.05
line.set_ydata(u_heat(x, t))
title.set_text(f't = {t:.2f}')
return line, title
ani = FuncAnimation(fig, animate, frames=120, interval=50, blit=True)
plt.show()Exercise 2.1.7. Numerically verify d'Alembert's formula: for \(c=1\), \(\phi(x)=\exp(-10x^2)\), \(\psi(x)=0\), plot the solution at \(t=0,0.5,1,2\) and confirm the two pulses propagate at unit speed.
Exercise 2.1.8. Compute and plot the Poisson kernel for the unit disk:
$$ P(r,\theta) = \frac{1-r^2}{1-2r\cos\theta+r^2}, $$for \(r=0.3,0.6,0.9,0.99\). Observe how the kernel concentrates as \(r\to 1\).
Exercise 2.1.9 (Agentic). Design a decision-tree agent that takes as input the coefficients \((A,B,C)\) of a second-order PDE plus a list of boundary/initial conditions, and outputs: (i) the classification, (ii) whether the problem is well-posed, and (iii) a recommended solution strategy. Test it on at least 10 distinct PDEs from your textbook.
Exercise 2.1.10 (Agentic). Implement a multi-armed bandit that selects the number of Fourier terms \(N\) to use in the heat-equation solution. The reward is \(-(\text{error} + \alpha\cdot N)\) where error is the \(L^2\) norm against a reference solution with \(N=500\) and \(\alpha\) is a computational cost penalty. Run 1000 episodes and plot the average \(N\) chosen over time.
Q1. The PDE \(u_{xx}+u_{yy}+u_{zz}=0\) is:
(a) Parabolic (b) Hyperbolic (c) Elliptic (d) None of the above
Q2. The discriminant \(\Delta=B^2-AC\) for the wave equation \(u_{tt}=c^2 u_{xx}\) (written as \(c^2 u_{xx}-u_{tt}=0\)) is:
(a) \(-c^2\) (b) \(0\) (c) \(c^2\) (d) \(1\)
Q3. Separation of variables applied to the heat equation on \([0,L]\) with Dirichlet BCs yields eigenvalues:
(a) \(\lambda_n = n\pi/L\) (b) \(\lambda_n = (n\pi/L)^2\) (c) \(\lambda_n = n^2\) (d) \(\lambda_n = e^{-n\pi/L}\)
Q4. In d'Alembert's formula, the domain of dependence for the point \((x_0,t_0)\) is:
(a) \([x_0-ct_0,\;x_0+ct_0]\) (b) \([0,\;ct_0]\) (c) All of \(\mathbb{R}\) (d) \(\{x_0\}\)
Q5. The Gibbs phenomenon states that near a jump discontinuity, the Fourier partial sums overshoot by approximately:
(a) 5% (b) 9% (c) 15% (d) The overshoot vanishes as \(N\to\infty\)
Q6. A function \(u\) harmonic in a bounded domain \(\Omega\) attains its maximum:
(a) At the center of \(\Omega\) (b) On the boundary \(\partial\Omega\) (c) At a saddle point (d) Nowhere
Q7. The energy \(E(t)=\frac{1}{2}\int_0^L u^2\,dx\) for the heat equation with homogeneous Dirichlet BCs satisfies:
(a) \(E'(t)>0\) (b) \(E'(t)=0\) (c) \(E'(t)\leq 0\) (d) \(E'(t)\) oscillates
Q8. Which of the following is NOT one of Hadamard's three conditions for well-posedness?
(a) Existence (b) Uniqueness (c) Continuous dependence on data (d) Finite propagation speed
Q9. The Fourier sine coefficient \(b_n\) of \(f(x)=1\) on \([0,\pi]\) is:
(a) \(\frac{2}{n\pi}[1-(-1)^n]\) (b) \(\frac{1}{n}\) (c) \(0\) (d) \(\frac{2}{\pi}\)
Q10. The Tricomi equation \(y\,u_{xx}+u_{yy}=0\) is elliptic when:
(a) \(y>0\) (b) \(y<0\) (c) \(y=0\) (d) Always
Q1: (c) Q2: (c) Q3: (b) Q4: (a) Q5: (b)
Q6: (b) Q7: (c) Q8: (d) Q9: (a) Q10: (a)
Objective: Build a Python program that accepts user-specified coefficients \((A,B,C)\), classifies the PDE, and — for the canonical types — solves and visualises the solution using separation of variables or d'Alembert's formula.
| Criterion | Excellent (90--100%) | Good (70--89%) | Needs Improvement (<70%) |
|---|---|---|---|
| Correctness | All three PDE types classified correctly; solutions match analytic benchmarks to \(10^{-6}\) relative error with sufficient Fourier terms. | Two of three types handled correctly; minor truncation issues. | Misclassification or incorrect solutions. |
| Efficiency | Vectorised Fourier sums; animation runs at 30+ fps; code profiled and bottlenecks documented. | Runs in reasonable time but uses explicit Python loops. | Unacceptably slow for \(N>20\) terms. |
| Robustness | Handles edge cases (zero coefficients, degenerate BCs) gracefully with informative error messages; input validation present. | Handles typical inputs; crashes on edge cases. | Crashes on valid inputs. |
| Reproducibility | Includes requirements.txt, docstrings, and a README with sample outputs. Random seeds fixed where applicable. |
Code runs but lacks documentation. | Missing dependencies or unclear usage. |