1. Steady heat conduction: Suppose T(x) is temperature and k(x) is thermal conductivity. Fourier's law is
   \[ q = -k(x)\,\nabla T . \]
   At steady state, with no internal heat generation,
   \[ \nabla \cdot q = 0 \qquad \Longrightarrow \qquad \nabla \dot k(x)\,\nabla T = 0 . \]
   This is the variable-coefficient Laplace equation. If k is constant, then
   \[ k\,\nabla^2 T = 0 \qquad \Longrightarrow \qquad \nabla^2 T = 0 . \]
2. Electrostatics in an inhomogeneous dielectric: Let ϕ(x) be the electric potential, measured in volts. Electrostatic potential (or electric potential) is the amount of work needed to move a unit positive test charge from a reference point (like infinity) to a specific point inside a static electric field:
   \[ \mathbf{E} = - \nabla \phi . \]
   Then the electric displacement field, denoted by 𝔃, is
   \[ \mathbf{D} = \epsilon(\mathbf{x})\,\mathbf{E} . \]
   Therefore,
   \[ \mathbf{D} = -\epsilon(\mathbf{x})\,\nabla \phi . \]
   We have Gauss's law:
   \[ \nabla \cdot \mathbf{D} = \rho , \]
   and, in a charge-free region (ρ = 0), and a material with spatially varying permittivity, ϵ = ϵ(x):
   \[ \nabla \cdot \mathbf{D} = 0 \qquad \Longrightarrow \qquad \nabla \dot \left( \epsilon (\mathbf{x}\,\nabla \phi \right) = 0 . \]
   Again, if ϵ is constant,
   \[ \nabla^2 \phi = 0 . \]
3. Steady diffusion Suppose c(x) is concentration and D(x) is a spatially varying diffusion coefficient. Fick's law gives
   \[ \mathbf{J} = - \mathbf{D}(\mathbf{x})\,\nabla c . \]
   At steady state, with no sources or sinks,
   \[ \nabla \dot mathbf{J} = 0\qquad \Longrightarrow \qquad \nabla \dot \left( \mathbf{D}(\mathbf{x})\,\nabla c \right) = 0 . \]
   For constant 𝔃,
   \[ \nabla^2 c = 0 . \]
   So heat, electrostatics, and diffusion all have essentially the same mathematical structure.

4. Potential flow: For an incompressible, irrotational fluid,
   \[ \nabla \cdot \mathbf{v} = 0 , \qquad \nabla \times \mathbf{v} = 0 . \]
   If the flow is irrotational, we introduce a velocity potential
   \[ \mathbf{v} = \nabla\,\phi . \]
   Then incompressibility gives
   \[ \nabla \cdot \nabla \,\phi = 0 \qquad \iff \qquad \nabla^2 \phi = 0 . \]
   This is the constant-coefficient Laplace equation.
5. Laplace's equation repeatedly appears whenever you have a conserved flux proportional to a gradient. The common structure is
   \[ \nabla \cdot \left( K(x) \,\nabla u \right) = 0 , \]
   where K may be a scalar k(x) or, more generally, a matrix/tensor.

   Conservation law Let's start with the most fundamental derivation. Suppose u(x,y) is some scalar field and its flux is

   \[ \mathbf{J} = -k\,\nabla u , \]
   where k is a constant. In two dimensions,
   \[ \mathbf{J} = -k\,\begin{pmatrix} u_x \\ u_y \end{pmatrix}\nabla u , \]
   Now impose steady-state conservation:
   \[ \nabla \cdot \mathbf{J} = 0 . \]
   Therefore,
   \[ \nabla \cdot \left( -k\,\nabla u \right) = 0 . \]
   Since k is constant,
   \[ -k\,\nabla \cdot \left( \nabla u \right) = 0 \qquad \Longrightarrow \qquad u_{xx} + u_{yy} = 0 . \]
   That's the two-dimensional Laplace equation.
6. Variable thermal conductivity Now let's make things more interesting. Let T(x,y) be temperature and suppose thermal conductivity varies in space: k = k(x,y). Fourier's law gives
   \[ q = \begin{pmatrix} -k\,T_x \\ -k\,T_y \end{pmatrix} . \]
   At steady state, assuming there is no heat generation,
   \[ \nabla \cdot q = 0 . \]
   Thus,
   \[ \frac{\partial}{\partial x} \left( -k\,T_x \right) + \frac{\partial}{\partial y} \left( -k\,T_y \right) = 0 \qquad\iff \qquad \frac{\partial}{\partial x} \left( k\,T_x \right) + \frac{\partial}{\partial y} \left( k\,T_y \right) = 0 \]
   In succinct notation,
   \[ \nabla \cdot \left( k\,\nabla T \right) = 0 . \]
   Note that the principal symbol for variable k and constant k is the same,
   \[ \sigma (L) = - k\, \| \xi \|^2 . \]

7. Variational principle: Suppose we define an energy functional
   \[ E[u] = \frac{1}{2}\,\int_{\Omega} \,k(\mathbf{x}) \left\vert \nabla u \right\vert^2 \,{\text d}V . \]
   We ask: Which function u minimizes this energy?

   Take a small variation

   \[ u \,\to\, u + \\varepsilon \,v , \]
   where v = 0 on the boundary. Then
   \[ E[u + \varepsilon\, v] = \frac{1}{2}\,\int_{\Omega} \,k(\mathbf{x}) \left\vert \nabla u + \varepsilon\, \nabla v\right\vert^2 \,{\text d}V \]
   Expand:
   \[ E[u + \varepsilon\, v] = \frac{1}{2}\,\int_{\Omega} \,k(\mathbf{x}) \left( |\nabla u|^2 + 2\varepsilon \,\nabla u \cdot \nabla v + \varepsilon^2 \, |\nabla v |^2 \right) \,{\text d}V . \]
   For an extremum, \( \displaystyle \quad \left. \frac{{\text d}E}{{\text d}\varepsilon} \right\vert_{\varepsilon = 0} = 0 . \quad \) Therefore,
   \[ \int_{\Omega}\, k\,\nabla u \cdot \nabla v \,{\text d}V = 0 . \]
   Now integrate by parts
   \[ \int_{\Omega}\, k\,\nabla u \cdot \nabla v \,{\text d}V = - \int_{\Omega}\, v\,\nabla \cdot \left( k\,\nabla u \right) \,{\text d}V + \int_{\partial\Omega}\, vk\,\nabla u \cdot \hat{\bf n} \,{\text d}S . \]
   Since v = 0 on the boundary,
   \[ \int_{\Omega}\, v\,\nabla \cdot \left( k\,\nabla u \right) \,{\text d}V = 0 . \]
   Since this must hold for arbitrary v,
   \[ \nabla \cdot \left( k\,\nabla u \right) = 0 . \]
   And for constant k,
   \[ \nabla^2 u = 0 . \]
   We have now obtained Laplace's equation without starting from a conservation law.

8. Anisotropic medium: Now for the most general and interesting case, suppose the flux does not necessarily point opposite to the gradient. Instead,
   \[ \mathbf{J} = - \mathbf{K}(\mathbf{x})\,\nabla u , \]
   where
   \[ \mathbf{K} = \begin{pmatrix} k_{xx} & k_{xy} \\ k_{yx} & k_{yy} \end{\pmatrix} . \]
   This happens, for example, in an anisotropic material where conductivity depends on direction. Then conservation gives
   \[ \nabla \cdot \mathbf{J} = 0 . \]
   Therefore,
   \[ \nabla \cdot \left( \mathbf{K} (\mathbf{x} )\,\nabla u \right) = 0 . \]
   In two dimensions,
   \[ \mathbf{K}\,\nabla u = \begin{pmatrix} k_{xx} u_x + k_{xy} u_y \\ k_{yx} u_x + k_{yy} u_y \end{pmatrix} . \]
   Thus,
   \[ \frac{\partial}{\partial x} \left( k_{xx} u_x + k_{xy} u_y \right) + \frac{\partial}{\partial y} \left( k_{yx} u_x + k_{yy} u_y \right) = 0 . \]
   Its principal symbol is
   \[ \sigma_2 (\xi ) = - \left( k_{xx} \xi_1^2 + k_{xy} \xi_1 \xi_2 + k_{yx} \xi_1 \xi_2 + k_{yy} \xi_2^2 . \]
   If all coefficients are constant, expanding gives
   \[ k_{xx} u_{xx} + \left( k_{xy} + k_{yx} \right) u_{xy} + k_{yy} u_{yy} = 0 . \]
   If the conductivity tensor is symmetric, k_xy = k_yx, then
   \[ k_{xx} u_{xx} + 2\,k_{xy} \, u_{xy} + k_{yy} u_{yy} = 0 . \]
   If K = kI, then
   \[ k_{xx} = k_{yy} , \qquad k_{xy} = 0 , \]
   so we recover
   \[ k\,u_{xx} + k\,u_{yy} = 0 \qquad \iff \qquad u_{xx} + u_{yy} = 0 . \]
      ■
   End of Example 1

The biharmonic operator Δ² has the principal symbol \[ \sigma_4 \left( \Delta^2 \right)\,( \xi ) = |\xi |^4 = \left( \xi_1^2 + \xi_2^2 + \cdots + \xi_n^2 \right)^2 . \] Since it is nonzero for ξ ≠ 0, the biharmonic operator is elliptic.

The biharmonic equation can be written in two dimensions as

Here are some applications where the appearance of the fourth-order operator is genuinely meaningful.

1. Thin elastic plates: For a thin Kirchhoff–Love plate, the transverse displacement u(x,y) satisfies
   \[ D\,\Delta^2 u = f , \]
   where D = Eh³/[12(1−ν²)] is the flexural rigidity. Here E is Young's modulus h is elastic thickness, ν is Poisson's ratio.

   The interesting part is not merely that a plate bends. The biharmonic operator comes from the fact that the elastic energy depends on second derivatives of displacement:

   \[ E[u] \,\sim\, \int_{\Omega} , \left[ (\Delta\, u)^2 + (1-\nu )\left( 2\,u_{xy}^2 - u_{xx} u_{yy} \right) \right] {\text d}x{\text d}y . \]
   Taking the Euler–Lagrange equation of this functional produces a fourth-order PDE.

   This becomes particularly nontrivial for free boundaries. You don't simply prescribe u and ∂n​u; the natural boundary conditions involve bending moments and Kirchhoff shear forces. Thus the biharmonic problem encodes both the deformation and the mechanics at the boundary.

2. Stokes flow and the stream function: A steady, incompressible, 2D Stokes flow (also known as creeping flow) describes the motion of fluids where viscous forces completely dominate advective inertial forces, occurring at very low Reynolds numbers (Re ≪ 1) where the nonlinear convective acceleration terms vanish. The equation of motion for Stokes flow (George Stokes, 1842--1850) can be obtained by linearizing the steady state Navier–Stokes equations:
   \[ - \nabla p + \mu\,\Delta u = 0 , \qquad \nabla \cdot u = 0 . \]
   introduce a stream function ψ:
   \[ u = \psi_y , \qquad v = - \psi_x , \]
   he incompressibility condition is then automatically satisfied. Taking the curl of the Stokes equations eliminates pressure and gives
   \[ \Delta \omega = 0 , \]
   where the vorticity is
   \[ \omega = v_x - u_y = - \Delta \psi . \]
   Consequently,
   \[ \Delta^2 \psi = 0 . \]
   This is important in problems such as flow past a cylinder, flow through narrow geometries, lubrication-type problems, and the Stokes paradox in two dimensions.
3. Airy stress function in 2D elasticity: Another beautiful appearance of the biharmonic operator is in linear elasticity. For a 2D elastic body with no body forces, equilibrium requires
   \[ \begin{split} \sigma_{xx,x} + \sigma_{xy,y} &= 0 , \\ \sigma_{xy,x} + \sigma_{yy,y} &= 0 , \end{split} \]
   where σ is the Cauchy stress tensor, ε is the infinitesimal strain tensor, u is the displacement vector. Introduce the Airy stress function ϕ:
   \[ \sigma_{xx} = \frac{\partial^2 \phi}{\partial y^2} , \qquad \sigma_{xy} = - \frac{\partial^2 \phi}{\partial x\,\partial y} \qquad \sigma_{yy} = \frac{\partial^2 \phi}{\partial x^2} . \tag{A.1} \]
   Equilibrium is now satisfied identically. Then compatibility of the strains, together with Hooke's law, gives
   \[ \Delta^2 \phi = \nabl;a^4 \phi = 0 \]
   for homogeneous isotropic elasticity without body forces. Any function φ satisfying Δ²u = 0 is called a bihar- monic function. A harmonic function f satisfies ∇²f = 0. Thus, if f is harmonic, it is also biharmonic. However, the converse is not true. Once the Airy stress function is known, the stresses can be obtained by Eq.(A.1).

   What's particularly interesting here is that the biharmonic function is not the displacement. It is a scalar potential for the stress field.

4. Hele–Shaw flow and viscous fingering: In a Hele–Shaw cell, a viscous fluid flows between two closely spaced plates. Under the depth-averaged approximation, the velocity is governed by Darcy's law,
   \[ u = - \frac{b^2}{12\,\mu}\,\nabla p . \]
   Darcy's law is a simple equation in physics and engineering that tells us how fast water or other liquids move through tiny spaces inside rocks, sand, or soil. French engineer Henri Darcy discovered it in 1856 by studying water flow through sand filter.

   Since the fluid is incompressible,

   \[ \nabla \cdot \mathbf{u} = 0 , \]
   so the pressure is harmonic:
   \[ \Delta p = 0 . \]
   That sounds like an ordinary Laplace problem. But introducing a stream function and considering the velocity field leads naturally to higher-order formulations; in particular, biharmonic structures appear in the boundary-value formulation of the velocity/stream function.

   This becomes especially interesting when studying Saffman–Taylor viscous fingering (Saffman–Taylor instability), where the interface between two fluids is unstable. The harmonic pressure field couples to a moving free boundary, so the resulting problem is much richer than simply solving Δu = 0.

5. Surface interpolation / thin-plate splines: This is an application from approximation theory and machine learning rather than classical physics. Suppose you have scattered measurements
   \[ \left( x_i , \, y_i\, z_i \right) \]
   and want a smooth surface u(x,y) passing through them. One can minimize the bending energy
   \[ J[u] = \,\iint_{\mathbb{R}^2} \,\left( u_{xx}^2 + 2\,u_{xy}^2 + u_{yy}^2 \right) {\text d}x {\text d}y \]
   subject to
   \[ u\left( x_i , y_i \right) = z_i . \]
   The Euler–Lagrange equation away from the interpolation points is
   \[ \Delta^2 u = 0 . \]
   The fundamental solution of the 2D biharmonic operator is, up to constants,
   \[ r^2 \ln r . \]
   Hence the famous thin-plate spline representation
   \[ u(x,y) = a_0 + a_1 x + a_2 y + \sum_i \lambda_i \| x - x_i \|^2 \ln \| x - x_i \| . \]
   This is a particularly nice example because the Green's function of the biharmonic operator directly becomes the interpolation kernel.
6. Gradient flows and phase-field models: There is a more modern connection to Cahn–Hilliard dynamics that describes the process of phase separation, spinodal decomposition, by which the two components of a binary fluid spontaneously separate and form domains pure in each component. For a free energy such as
   \[ F[c] = \int_{\Omega} \left( f(c) + \frac{\kappa}{2}\,\left\vert \nabla c \right\vert^2 \right) {\text d}x , \]
   where c is the concentration of the fluid, the chemical potential is
   \[ \mu = f' (c) - \kappa\,\Delta c . \]
   Conserved dynamics takes the form
   \[ c_t = M\,\Delta \mu . \]
   Therefore,
   \[ c_t = M\,\Delta f'(c) - M\,\kappa \,\Delta^2 c . \]
   So the biharmonic operator appears because you have a conserved gradient flow whose chemical potential already contains a Laplacian.

   For constant mobility and a quadratic f, this essentially give

   \[ c_t = - M\,\kappa \,\Delta^2 c + \cdots . \]
   This is important for spinodal decomposition, phase separation, coarsening, and interfacial dynamics.

7. Plate/shell theory and higher-gradient mechanics: A more sophisticated generalization occurs in strain-gradient and couple-stress elasticity. Classical elasticity usually involves first derivatives of displacement in the strain energy. But if the material energy depends on gradients of strain,
   \[ E[u] \,\sim\,\int_{\Omega} \left( \left\vert \nabla u \right\vert^2 + \ell^2 \left\vert \nabla^2 u \right\vert^2 \right) {\text d}x , \]
   then variation produces operators such as
   \[ - \Delta u + \ell^2 \Delta^2 u . \]
   Thus, biharmonic terms arise naturally when the material possesses an intrinsic length scale ℓ.

Let u : ℝ² → ℝ² be a plane displacement field. The linear elasticity operator is actually the Lamé operator In homogeneous and isotropic materials \[ L(\mathbf{u}) = \mu\,\Delta\,\mathbf{u} + \left( \lambda + \mu \right) \nabla \left( \nabla \cdot \mathbf{u} \right) , \] where λ and μ are positive Lamé parameters, two material-dependent quantities. It is written in vector form, so the three dimensional case is treated similarly. The principal symbol is a 3 × 3 matrix: \[ \sigma_2 (L)(\mathbf{x}, \xi ) = \mu \left\vert \xi \right\vert^2 \mathbf{I} + \left( \mu + \lambda \right) \xi \otimes \xi , \] where

• I is the identity matrix,
• ξ ⊗ ξ is the rank one matrix with entries ξ_iξ_j.

Let us show ellipticity of Lamé's operator. We work with ℝⁿ for n = 2, 3. We want to determine eigenvalue of matrix \( \displaystyle \quad \mathbf{S}(\xi ) = \sigma_2 (L) (\xi ) = \mu \|\xi \|^2 + \left( \lambda + \mu \right) \xi\otimes \xi . \quad \)

For fixed ξ ≠ 0, we consider

• the longitudinal subspace \[ V_{\parallel} = \left\{ \mathbf{v} \in \mathbb{R}^n : \quad \mathbf{v} = \alpha\xi \mbox{ for some } \alpha \in\mathbb{R} \right\} , \]
• the transverse subspace \[ V_{\perp} = \left\{ \mathbf{w} \in \mathbb{R}^n \ : \quad \mathhbf{w} \cdot \xi = 0 \mbox{ for } \xi \in V_{\parallel} \right\} . \]
<

Take w ∈ V_∥, so w · ξ = 0. WE compute

Take v ∈ V_∥, so v = α ξ for some α ∈ ℝ. Compute

Hence, the eigenvalues of Lamé's symbol are

• μ|ξ|² with multiplicity n − 1,
• (λ + 2μ) |ξ|² with multiplicity 1.

The traction vector on the boundary r = R is

We consider a thin plate occupying a circular domain whose boundary is given by the cylindrical surface r = R. The plate is modeled according to Kirchhoff–Love plate theory, in which the deformation of the three-dimensional plate is described in terms of the displacement of its two-dimensional mid-surface. In particular, the transverse displacement is represented by a scalar function w(r,θ) defined on the circular domain. The governing equation for w is a fourth-order biharmonic equation, with boundary conditions prescribed at r = R.

For an isotropic plate of flexural rigidity

The shear forces in terms of deflection are

There is another mathematically elegant possibility of boundary conditions

References:

• Gazzola, F., Grunau, H.-C., Sweers, G., Polyharmonic Boundary Value Problems, Springer, 2010, doi: 10.1007/978-3-642-12245-3

  discusses different boundary-value problems for the biharmonic and more general polyharmonic operators.

• Timoshenko, S., Woinowsky-Krieger, S. Theory of Plates and Shells, McGraw-Hill College, 1959. 2nd ed.,

  which treats the classical plate boundary conditions and circular plates.

• Plates from 3D elasticity,

  for a more modern treatmen.

The traction operator is

1. Krainer, T., (2005) Elliptic boundary problems on manifolds with polycylindrical ends, arXiv:math/0508516 [math.AP]
2. Krupchyk, K. & Tuomela, J., (2006) The Shapiro–Lopatinskij Condition for Elliptic Boundary Value Problems, London Mathematical Society
3. Lopatinskii, Ya.B., (1953) On a method of reducing boundary problems for a system of differential equations of elliptic type to regular integral equations Journal: Ukraine. Mat. Zh., 5 (1953), 123–151.
5. Shapiro, Z.Ya,, (1953) On general boundary problems for equations of elliptic type, Journal: Izvestiya Akad. Nauk SSSR. Ser. Mat., 17 (1953), 539–562.