Solving a partial differential equation (PDE) means finding a function that satisfies the equation along with its boundary and initial conditions. There is no single method that works for all PDEs. The right approach depends on the type of equation, the geometry of the problem, and whether you need an exact formula or a numerical approximation. This guide walks through the major techniques, starting with how to classify what you’re dealing with.
Classify the Equation First
Before choosing a solution method, you need to know what type of PDE you have. For second-order linear PDEs, classification comes down to a simple discriminant test. The general form is Auxx + Buxy + Cuyy + lower-order terms = 0, and the discriminant is B² − 4AC. If the discriminant is positive, the equation is hyperbolic. If it equals zero, it’s parabolic. If it’s negative, it’s elliptic.
Each type corresponds to fundamentally different physical behavior. Hyperbolic equations describe wave propagation, with information traveling along specific paths called characteristic curves. Parabolic equations describe diffusion processes like heat spreading through a material, with a single family of characteristics. Elliptic equations describe steady-state situations like temperature distributions that have settled into equilibrium, and they have no characteristic curves at all. The three canonical examples are:
- Wave equation (hyperbolic): utt = a²∇²u
- Heat equation (parabolic): ut = k∇²u
- Laplace’s equation (elliptic): ∇²u = 0
This classification tells you which solution methods will work. Hyperbolic equations respond well to the method of characteristics. Parabolic and elliptic equations on bounded domains often yield to separation of variables. Problems on infinite domains call for integral transforms.
Boundary and Initial Conditions
A PDE alone doesn’t have a unique solution. You need additional conditions that pin down one specific answer. Initial conditions tell you the state at time zero (the starting temperature profile, the initial shape of a vibrating string). Boundary conditions tell you what happens at the edges of your domain for all time.
The three standard types of boundary conditions are Dirichlet, Neumann, and Robin. Dirichlet conditions fix the value of the solution at the boundary, like holding the ends of a rod at specific temperatures: u(0, t) = 0. Neumann conditions fix the derivative at the boundary, which in heat problems corresponds to controlling the heat flux: ux(0, t) = 0 means the boundary is insulated. Robin conditions combine both, linking the value and its derivative at the boundary: ux(L, t) + κu(L, t) = 0. This models situations where heat escapes at a rate proportional to the temperature difference with the surroundings.
For a PDE problem to be well-posed in the classical sense (often called Hadamard well-posedness), three things must hold: a solution exists, it is unique, and it depends continuously on the data. That last condition means small changes in your boundary or initial conditions produce only small changes in the solution, not wild swings. If any of these fail, you likely have the wrong type or number of conditions for your equation.
Separation of Variables
Separation of variables is the most widely taught analytical method and works on linear PDEs with homogeneous boundary conditions on regular geometric domains (rectangles, circles, spheres). The core idea is to guess that the solution can be written as a product of functions, each depending on only one variable. For a function of x and t, you assume u(x, t) = X(x)T(t).
The process has four steps. First, substitute the product form into the PDE. Because each side of the resulting equation depends on a different variable, both sides must equal the same constant (the separation constant λ). This splits your PDE into two ordinary differential equations, which are far easier to solve.
Second, apply the boundary conditions. These restrict which values of λ are allowed, typically producing a discrete indexed set λ₁, λ₂, λ₃, and so on. Each allowed λ gives you a “modal solution” that satisfies both the PDE and the boundary conditions.
Third, use superposition. Since the equation is linear, any sum of modal solutions is also a solution. The general solution is an infinite series combining all the modes, each multiplied by an unknown coefficient.
Fourth, determine the coefficients using the initial conditions. This usually involves matching the series to the initial data and using orthogonality of the basis functions (sines, cosines, or Bessel functions, depending on the geometry) to extract each coefficient individually. In practice, this means computing Fourier coefficients.
For Laplace’s equation on a rectangle, the same logic applies but with two spatial variables instead of one spatial and one time variable. The separated solutions take forms like sinh(nπx/b)sin(nπy/b), and the general solution is a series with coefficients chosen to match the boundary data.
Method of Characteristics
For first-order PDEs, the method of characteristics converts the PDE into a system of ordinary differential equations (ODEs) that you can solve with standard techniques. The idea is geometric: information in a first-order PDE propagates along specific curves in the x-t plane called characteristics. Along each curve, the PDE reduces to an ODE.
For a quasilinear first-order PDE of the form ut + c(x, t, u)ux = f(x, t, u) with initial data u(x, 0) = u₀(x), the recipe is straightforward. Pick a starting point x₀ on the initial line. Then solve the coupled ODE system: X'(t) = c(X, t, v) with X(0) = x₀, and v'(t) = f(X, t, v) with v(0) = u₀(x₀). The function X(t) traces out the characteristic curve, and v(t) gives the value of the solution along that curve.
Repeating this for every starting point x₀ sweeps out the full solution. The method works cleanly until characteristics cross, which can happen in nonlinear problems and signals the formation of a shock wave. At that point, you need additional theory (weak solutions, entropy conditions) to continue.
Integral Transforms
When your domain is infinite or semi-infinite, separation of variables breaks down because there are no boundary conditions to discretize the spectrum. Integral transforms handle this by converting the PDE in one variable into an ODE in the remaining variable.
The Fourier transform is the standard tool for problems on the entire real line. You transform with respect to the spatial variable x, which converts every x-derivative into multiplication by the transform variable ξ. The heat equation ut = kuxx on the whole line becomes ût = −kξ²û, a simple ODE in t with solution û(ξ, t) = ĝ(ξ)e−kξ²t, where ĝ is the transform of the initial data. Inverting the transform gives you the solution as a convolution with the heat kernel.
The same approach handles the wave equation, where the transformed solution involves sines and cosines in time: û(ξ, t) = ĝ(ξ)cos(c|ξ|t) + ĥ(ξ)sin(c|ξ|t)/(c|ξ|). Laplace’s equation on a half-plane transforms into an ODE in y with exponential solutions, and you discard the growing exponential to keep the solution bounded.
The Laplace transform (with respect to time rather than space) is particularly useful for problems with discontinuous forcing or complicated time-dependent boundary conditions. Multidimensional problems work the same way: you transform with respect to all spatial variables simultaneously, solve the resulting ODE, and invert.
Finite Difference Methods
When the geometry is complicated, the coefficients vary in space, or the equation is nonlinear, analytical methods often fail. Numerical methods approximate the solution on a grid of discrete points. Finite differences are the most intuitive approach: replace each derivative with an algebraic formula involving nearby grid values.
The first derivative can be approximated three ways. A forward difference uses the point ahead: (u(x+h) − u(x))/h. A backward difference uses the point behind: (u(x) − u(x−h))/h. Both are first-order accurate, meaning the error shrinks proportionally to the grid spacing h. A centered difference uses both neighbors: (u(x+h) − u(x−h))/(2h), which is second-order accurate, so halving the grid spacing cuts the error by a factor of four. The second derivative uses three points: (u(x−h) − 2u(x) + u(x+h))/h², also second-order accurate.
For the heat equation, the simplest scheme (forward in time, centered in space) computes each new time value explicitly from the current values: Uin+1 = Uin + (k/h²)(Ui−1n − 2Uin + Ui+1n). This is easy to code but only stable when the time step k is small enough relative to h² (specifically, k/h² ≤ 1/2). The Crank-Nicolson scheme averages the spatial derivatives between the current and next time step, giving second-order accuracy in both space and time with unconditional stability. The tradeoff is that you must solve a system of linear equations at each time step rather than just evaluating a formula.
For advection equations, stability depends on which direction information flows. Upwind methods use a backward difference in the direction of propagation, which naturally respects the physics and stays stable. Centered schemes for advection are unconditionally unstable unless modified (the Lax-Friedrichs method adds numerical diffusion to fix this).
Finite Element Methods
Finite element methods handle complex geometries and irregular domains more naturally than finite differences. The core idea is to reformulate the PDE in a “weak form” that relaxes the smoothness requirements on the solution.
The conversion happens in two steps. First, multiply the PDE by an arbitrary test function and integrate over the domain. Second, use integration by parts (Green’s formula) to shift one derivative from the solution onto the test function. For the Poisson equation −Δu = f, this transforms the requirement that the second derivative equals f everywhere into the requirement that ∫∇u·∇w dx = ∫fw dx for every test function w. This is the weak formulation.
The practical payoff is that you only need first derivatives of u, not second derivatives. You can then approximate u as a combination of simple basis functions (typically piecewise polynomials defined on triangles or tetrahedra that tile the domain). Substituting this approximation into the weak form and choosing each basis function as a test function converts the PDE into a system of linear algebraic equations, which you solve with standard matrix methods. Refining the mesh (using more, smaller elements) improves accuracy at the cost of a larger system to solve.
Software for Solving PDEs
For problems beyond textbook exercises, software handles the heavy computation. FEniCS is the most widely used open-source finite element library, available for both Python and Julia. You write the weak form of your problem in near-mathematical notation, and FEniCS handles meshing, assembly, and solving. It’s particularly strong for elliptic and parabolic problems on complex domains.
For finite difference approaches, Python’s SciPy provides basic ODE solvers that work after you discretize the spatial dimensions (the “method of lines”). Julia’s SciML ecosystem offers more specialized tools: MethodOfLines.jl automates the spatial discretization, while DiffEqOperators.jl gives finer control over finite difference schemes. Trixi.jl handles hyperbolic conservation laws with adaptive high-order methods, which is useful for fluid dynamics problems with shocks.
For spectral methods (which use global basis functions instead of local ones and converge extremely fast for smooth solutions), ApproxFun.jl and Python’s Dedalus project are strong choices. Physics-informed neural networks, which train a neural network to satisfy the PDE, are available through NeuralPDE.jl and PyTorch-based libraries, though they remain less reliable than traditional methods for most standard problems.
Choosing the Right Method
Linear PDE on a simple domain (rectangle, disk, sphere) with constant coefficients: try separation of variables first. You’ll get an exact series solution.
First-order PDE: use the method of characteristics. It gives exact solutions and reveals the geometry of how information propagates.
Linear PDE on an infinite domain: reach for Fourier or Laplace transforms. The transform converts the PDE to an ODE, and the main challenge shifts to computing the inverse transform.
Nonlinear PDE, variable coefficients, or irregular geometry: go numerical. Finite differences are simplest to implement and work well on rectangular grids. Finite elements handle complex shapes and are the standard in engineering simulation. For either approach, start with a coarse grid to verify your setup, then refine until the solution stops changing meaningfully.