Poisson Equation, Fundamental Solution, and Newtonian Potential

Intuition Beginner

The Laplace equation describes equilibrium with no sources: the temperature in a metal sheet whose interior is undisturbed, the gravitational potential in empty space, the electrostatic potential where no charge sits. The Poisson equation is the same equation with sources allowed. Wherever a mass density or charge density is present, the potential is no longer harmonic; it has a Laplacian equal to minus that density (up to constants). The Poisson equation is written in shorthand as Lu = -f, where L is the Laplacian and f is the source.

Newton's gravity is the original Poisson equation. The gravitational potential of a mass distribution obeys an equation of the form Lu = c times the mass density, where c is a positive physical constant (in the modern convention, the constant is 4 pi G with G the gravitational constant). Outside the mass, the potential is harmonic and the equation reduces to the Laplace equation. Inside the mass, the potential bends toward the matter, and the Laplacian records exactly how much.

Electrostatics works the same way. The electric potential V satisfies the Poisson equation with the source f given by minus the charge density divided by the permittivity of vacuum. In a charge-free region V is harmonic, but inside a region of accumulated charge the potential curves and the Laplacian picks up the charge density. The pattern repeats across physics: any time a field comes from a density, the field's potential obeys a Poisson equation.

The remarkable fact is that one explicit function does almost all the work. In three-dimensional space, the function 1 over the distance from the origin is the potential of a point mass placed at the origin. This single function is called the fundamental solution. Once we know it, we can build the potential of any continuous mass distribution by adding up the contributions of every small piece of mass, weighted by the fundamental solution centred at that piece's location. The construction is called the Newtonian potential, and it gives the answer to the Poisson equation directly as a single integral.

In symbols: the potential u at a point x is the total of the source f times the fundamental solution centred at x, summed over all source points y. In two dimensions the fundamental solution is the negative logarithm of the distance divided by 2 pi; in three or more dimensions it is a power of the distance with a dimensional constant out front. The dimensions matter because the fundamental solution captures the geometry of how potentials spread out from a point source, and that geometry changes with dimension.

The one-sentence takeaway: the Poisson equation is the Laplace equation with sources, the fundamental solution is the potential of a single point source, and the Newtonian potential builds the solution to the Poisson equation by adding up point-source contributions across the whole source distribution.

Visual Beginner

Picture a single grain of mass at the origin of three-dimensional space. The gravitational potential it produces is a function that decays as 1 over the distance: large negative values near the grain, fading toward zero far away. Now replace the single grain with a continuous cloud of mass spread through a region.

The potential at any point outside the cloud is the sum of the contributions of every grain in the cloud, each contributing its own 1 over the distance term centred at its own location. Inside the cloud, the same formula still holds, with the singularity tamed because the source has finite density and the integral converges. The picture below shows a point-source potential next to a smooth-source potential, with the Poisson equation written above each as the local statement that the Laplacian of the potential equals the local source density (up to sign).

The visual point is the bridge between the local picture (the Laplacian at one point equals the source there) and the global picture (the potential at one point equals an integral over all sources). The fundamental solution is the kernel that converts the local picture into the global picture.

Worked example Beginner

We compute the gravitational potential of a uniform-density solid ball of radius R and total mass M sitting at the origin in three-dimensional space, using the Newtonian-potential idea. The answer is famous: outside the ball the potential is the same as if all the mass were concentrated at the origin, and inside the ball the potential grows quadratically with distance from the centre.

Step 1. The density is ${\rho }_0=3M/(4\pi R^3)$ inside the ball and zero outside. The Poisson equation reads $\Delta u=-4\pi G{\rho }_0$ inside the ball and $\Delta u=0$ outside.

Step 2. By symmetry, the potential depends only on the radial distance $r$ from the origin. Write $u=u(r)$. The radial Laplacian in three dimensions is $u^{\prime \prime }+(2/r)u^{\prime }$.

Step 3. Outside the ball, solve $u^{\prime \prime }+(2/r)u^{\prime }=0$. The general radial harmonic function is $u(r)=A/r+B$. Requiring $u\to 0$ at infinity gives $B=0$.

Step 4. Inside the ball, solve $u^{\prime \prime }+(2/r)u^{\prime }=-4\pi G{\rho }_0$. A particular solution is $u_p(r)=-\frac{2\pi G{\rho }_0}{3}{r}^2$ (verify: $u_p^{\prime \prime }=-4\pi G{\rho }_0/3$ and $(2/r)u_p^{\prime }=-8\pi G{\rho }_0/3$, summing to $-4\pi G{\rho }_0$). Add the radial harmonic correction $C+D/r$; finite potential at the centre forces $D=0$.

Step 5. Match values and slopes at $r=R$. Continuity gives $A/R=C-\frac{2\pi G{\rho }_0}{3}{R}^2$. Continuous derivative gives $-A/R^2=-\frac{4\pi G{\rho }_0}{3}R$. Solving: $A=\frac{4\pi G{\rho }_0{R}^3}{3}=GM$ and $C=-\frac{2\pi G{\rho }_0{R}^2}{1}=-\frac{3GM}{2R}$, after simplification.

Step 6. Write the final answer. Outside ($r\ge R$): $u(r)=GM/r$. Inside ($r\le R$): $u(r)=\frac{GM}{2R^3}{r}^2-\frac{3GM}{2R}$.

What this tells us: outside a spherically symmetric mass, the gravitational potential is identical to that of a point mass at the centre with the same total mass. This is Newton's shell theorem, recovered as a direct consequence of the Poisson equation and the fundamental solution.

Check your understanding Beginner

Formal definition Intermediate+

Let $\Omega \subseteq {\mathbb{R}}^n$ be an open set and let $f:\Omega \to \mathbb{R}$ be a measurable function. The Poisson equation is the second-order linear PDE $$-\Delta u=f\text{on}\Omega ,$$ where $\Delta u={\sum }_{i=1}^n{\partial }^{2}u/\partial x_i^2$ is the Laplacian. A classical solution is a function $u\in C^2(\Omega )$ satisfying the equation pointwise; a distributional solution is a distribution $u\in {\mathcal{D}}^{\prime }(\Omega )$ such that $-⟨u,\Delta \phi ⟩=⟨f,\phi ⟩$ for every test function $\phi \in C_c^{\infty }(\Omega )$ ^{[Evans 2010 §2.2.1]}.

Definition (fundamental solution). The fundamental solution of $-\Delta$ on ${\mathbb{R}}^n$ is the function $$\Phi (x)=\{\begin{array}{ll}-\frac{1}{2\pi }\log |x|,& n=2,\\ \frac{1}{n(n-2){\omega }_n}|x{|}^{2-n},& n\ge 3,\end{array}$$ where ${\omega }_n$ is the volume of the unit ball in ${\mathbb{R}}^n$. Equivalently, $n{\omega }_n=|S^{n-1}|$ is the surface area of the unit sphere in ${\mathbb{R}}^n$, and the $n\ge 3$ formula can be rewritten as $\Phi (x)=1/((n-2)|S^{n-1}||x{|}^{n-2})$.

The fundamental solution satisfies $\Delta \Phi =0$ on ${\mathbb{R}}^n\setminus \{0\}$ and $-\Delta \Phi ={\delta }_0$ in the sense of distributions on ${\mathbb{R}}^n$, where ${\delta }_0$ is the Dirac mass at the origin. The constant is fixed by requiring the surface-flux normalisation $-{\int }_{\partial B_r(0)}{\partial }_{\nu }\Phi dS=1$ for every $r>0$.

Definition (Newtonian potential). Given $f\in C_c^2({\mathbb{R}}^n)$ (compactly supported $C^2$), the Newtonian potential of $f$ is the convolution $$u(x)=(\Phi \ast f)(x)={\int }_{{\mathbb{R}}^n}\Phi (x-y)f(y)dy.$$ The construction extends to $f\in L^1({\mathbb{R}}^n)$ with compact support, and via further density arguments to general $L^p$ data with appropriate decay.

Spherical-symmetric reading. The Newtonian potential at $x$ is the integral over source space of the source density at $y$ multiplied by the fundamental-solution kernel $\Phi (x-y)$, which depends only on the distance $|x-y|$. The kernel captures how a unit point source at $y$ contributes to the potential at $x$; the integral superposes these contributions linearly. The linearity is essential: the Poisson equation is linear, so the response to a sum of sources is the sum of the responses.

Counterexamples to common slips Intermediate+

• Convolution with the fundamental solution does not always give a classical solution. For $f$ merely continuous, the Newtonian potential $u=\Phi \ast f$ is $C^1$ but not necessarily $C^2$; the second derivatives can fail to exist pointwise even though the equation $-\Delta u=f$ holds distributionally. The $C^2$ hypothesis on $f$ in Theorem (Newtonian potential) below is the cleanest sufficient condition; weaker hypotheses require Hölder-continuity (Schauder theory) or $L^p$-theory (Calderón-Zygmund).

• The Newtonian potential is not the unique solution. On ${\mathbb{R}}^n$ for $n\ge 3$, the Newtonian potential is the unique solution decaying at infinity, but adding any harmonic function gives another solution. On ${\mathbb{R}}^2$, the Newtonian potential has logarithmic growth at infinity and uniqueness requires specifying the growth rate; a free additive constant is unavoidable.

• The fundamental solution in two dimensions has a sign change. The two-dimensional fundamental solution $-\log |x|/(2\pi )$ is negative for $|x|>1$ and positive for $|x|<1$. In contrast, in dimension $n\ge 3$ the fundamental solution $|x{|}^{2-n}/(n(n-2){\omega }_n)$ is positive everywhere outside the origin. Sign-tracking mistakes when computing potentials in different dimensions are common.

• Distributional differentiation does not commute with pointwise inequalities. The function $|x|$ has classical derivative $\mathrm{sgn}(x)$ away from zero, and distributional derivative the same thing, but the second distributional derivative is $2{\delta }_0$, not zero. The same care is needed for the fundamental solution: $\Delta \Phi =0$ pointwise on ${\mathbb{R}}^n\setminus \{0\}$ and $-\Delta \Phi ={\delta }_0$ distributionally; the two statements are compatible because the distributional equation picks up the singularity at the origin that the pointwise equation cannot see.

Key theorem with proof Intermediate+

Theorem (Newtonian potential solves the Poisson equation; Poisson 1813). Let $f\in C_c^2({\mathbb{R}}^n)$. Define $u:{\mathbb{R}}^n\to \mathbb{R}$ by $$u(x)={\int }_{{\mathbb{R}}^n}\Phi (x-y)f(y)dy.$$ Then $u\in C^2({\mathbb{R}}^n)$ and $-\Delta u=f$ on ${\mathbb{R}}^n$ ^{[Poisson 1813]}.

Proof. Change variables in the integral via the substitution $z=x-y$, giving $$u(x)={\int }_{{\mathbb{R}}^n}\Phi (z)f(x-z)dz.$$ This form moves the $x$-dependence onto $f$, which is $C^2$ with compact support, while $\Phi$ is locally integrable. We exploit this to justify differentiation under the integral sign.

Step 1 (first derivatives). For each $i\in \{1,\dots ,n\}$ and each $x\in {\mathbb{R}}^n$: $$\frac{u(x+h{e}_i)-u(x)}{h}={\int }_{{\mathbb{R}}^n}\Phi (z)\frac{f(x+h{e}_i-z)-f(x-z)}{h}dz.$$ Since $f\in C_c^1$, the difference quotient $h^{-1}[f(x+h{e}_i-z)-f(x-z)]$ converges uniformly in $z$ to $\partial f/\partial x_i(x-z)$ as $h\to 0$, and is uniformly bounded by $\Vert \nabla f{\Vert }_{\infty }$. The bound is supported on a fixed compact set independent of $h$ (the support of $f$ shifted by $x$ together with a ball of radius $|h|$, which stays in a compact set). The dominated convergence theorem 02.07.04 applies, giving $${\partial }_{i}u(x)={\int }_{{\mathbb{R}}^n}\Phi (z){\partial }_{i}f(x-z)dz.$$ Repeating gives ${\partial }_i{\partial }_{j}u(x)=\int \Phi (z){\partial }_i{\partial }_{j}f(x-z)dz$ for every $i,j$, with the same justification using $f\in C_c^2$. So $u\in C^2({\mathbb{R}}^n)$.

Step 2 (apply the Laplacian). Summing over $i$: $$\Delta u(x)={\int }_{{\mathbb{R}}^n}\Phi (z){\Delta }_{x}f(x-z)dz={\int }_{{\mathbb{R}}^n}\Phi (z){\Delta }_{y}f(x-z)dz,$$ where the last equality uses that ${\Delta }_{x}f(x-z)={\Delta }_{y}f(x-z){|}_{y=x-z}$ since the Laplacian is invariant under the sign of the variable.

Step 3 (isolate the singularity). The integrand $\Phi (z){\Delta }_{y}f(x-z)$ is singular at $z=0$ where $\Phi$ blows up. Split the integration domain into $B_{\epsilon }=B_{\epsilon }(0)$ and ${\mathbb{R}}^n\setminus B_{\epsilon }$ for some small $\epsilon >0$: $$\Delta u(x)=\underset{I_{\epsilon }(x)}{\underbrace{{\int }_{B_{\epsilon }}\Phi (z)\Delta f(x-z)dz}}+\underset{J_{\epsilon }(x)}{\underbrace{{\int }_{{\mathbb{R}}^n\setminus B_{\epsilon }}\Phi (z)\Delta f(x-z)dz}}.$$

For the inner piece: $|\Phi (z)|$ is bounded by $C\log (1/|z|)$ in $n=2$ and by $C|z{|}^{2-n}$ in $n\ge 3$, both locally integrable. So $|I_{\epsilon }(x)|\le \Vert \Delta f{\Vert }_{\infty }{\int }_{B_{\epsilon }}|\Phi (z)|dz\to 0$ as $\epsilon \to 0$.

Step 4 (integrate by parts in the outer piece). On ${\mathbb{R}}^n\setminus B_{\epsilon }$, the function $\Phi$ is smooth. Apply Green's identity in the form ${\int }_D(\Phi \Delta f-f\Delta \Phi )={\int }_{\partial D}(\Phi {\partial }_{\nu }f-f{\partial }_{\nu }\Phi )dS$ on $D={\mathbb{R}}^n\setminus B_{\epsilon }$ (with the integration variable $z$ and $f$ shorthand for $f(x-z)$; the chain rule gives ${\Delta }_z[f(x-z)]=(\Delta f)(x-z)$). Since $\Delta \Phi =0$ on $D$ and $f(x-z)$ has compact support in $z$ (so the boundary contribution at infinity vanishes), the only surviving boundary terms are on $\partial B_{\epsilon }$: $$J_{\epsilon }(x)={\int }_{\partial B_{\epsilon }(0)}\left(\Phi (z){\partial }_{\nu }f(x-z)-f(x-z){\partial }_{\nu }\Phi (z)\right)dS(z),$$ where $\nu$ is the unit normal to $\partial B_{\epsilon }$ pointing into $B_{\epsilon }$ (outward from $D$).

Step 5 (boundary integrals as $\epsilon \to 0$). On $\partial B_{\epsilon }(0)$, $|\Phi |$ is comparable to ${\epsilon }^{2-n}$ (logarithmic in $n=2$) and the surface area is $|S^{n-1}|{\epsilon }^{n-1}$. So $$\left|{\int }_{\partial B_{\epsilon }}\Phi (z){\partial }_{\nu }f(x-z)dS\right|\le \Vert \nabla f{\Vert }_{\infty }\cdot |\Phi {|}_{\partial B_{\epsilon }}\cdot |S^{n-1}|{\epsilon }^{n-1}\to 0$$ as $\epsilon \to 0$ (the product ${\epsilon }^{2-n}\cdot {\epsilon }^{n-1}=\epsilon$ in $n\ge 3$; logarithmic in $n=2$ times $\epsilon$).

For the second boundary term, compute ${\partial }_{\nu }\Phi$ on $\partial B_{\epsilon }$. The inward normal is $-z/|z|=-z/\epsilon$. In $n\ge 3$: $$\Phi (z)=\frac{|z{|}^{2-n}}{n(n-2){\omega }_n},\nabla \Phi (z)=-\frac{|z{|}^{-n}z}{n{\omega }_n},{\partial }_{\nu }\Phi (z)=-\nabla \Phi (z)\cdot \frac{z}{\epsilon }=\frac{1}{n{\omega }_n{\epsilon }^{n-1}}.$$ In $n=2$, the same computation yields ${\partial }_{\nu }\Phi (z)=1/(2\pi \epsilon )$. In either case ${\partial }_{\nu }\Phi =1/(|S^{n-1}|{\epsilon }^{n-1})$.

So $${\int }_{\partial B_{\epsilon }(0)}f(x-z){\partial }_{\nu }\Phi (z)dS(z)=\frac{1}{|S^{n-1}|{\epsilon }^{n-1}}{\int }_{\partial B_{\epsilon }(0)}f(x-z)dS(z)\to f(x)$$ as $\epsilon \to 0$ by continuity of $f$ at $x$.

Step 6 (combine). Putting the pieces together: $$\Delta u(x)=\underset{\epsilon \to 0}{\lim }(I_{\epsilon }(x)+J_{\epsilon }(x))=0+(0-f(x))=-f(x).$$ So $-\Delta u=f$ on ${\mathbb{R}}^n$.

$\square$

Bridge. The Newtonian-potential construction builds toward 02.13.03 the heat equation, where the same convolution recipe applied to the heat-kernel fundamental solution gives the solution of the inhomogeneous Cauchy problem, and appears again in 02.13.01 the Laplace equation through the Green-function representation on bounded domains. The central insight is that the singular kernel $\Phi$ exactly compensates the boundary-of-a-shrinking-ball contribution in Green's identity, and this is exactly the foundational reason every constant-coefficient linear PDE on ${\mathbb{R}}^n$ admits an explicit solution operator: convolution with the fundamental solution. The bridge is between the local pointwise PDE $-\Delta u=f$ and the global integral-representation $u=\Phi \ast f$, identifying the differential operator $-\Delta$ on ${\mathbb{R}}^n$ with the inverse-by-convolution map on a suitable function space. Putting these together, the Riesz potential framework ^{[Riesz 1923]} generalises the construction to fractional powers $(-\Delta {)}^{-\alpha /2}$, and the Hörmander-Ehrenpreis-Malgrange existence theorem extends it to every non-zero constant-coefficient operator.

Exercises Intermediate+

Advanced results Master

The modern theory of the Poisson equation organises around five pillars: the fundamental-solution apparatus and Newtonian potential, the Green-function and layer-potential formulation of the Dirichlet and Neumann problems, the Riesz-potential and fractional-Laplacian framework, the Calderón-Zygmund regularity theory, and the general-existence apparatus for constant-coefficient operators (Hörmander-Ehrenpreis-Malgrange).

Theorem 1 (Riesz potentials; M. Riesz 1923). For $0<\alpha <n$, the Riesz potential of order $\alpha$ is the convolution operator $$I_{\alpha }f(x)=\frac{1}{{\gamma }_n(\alpha )}{\int }_{{\mathbb{R}}^n}\frac{f(y)}{|x-y{|}^{n-\alpha }}dy,{\gamma }_n(\alpha )=\frac{2^{\alpha }{\pi }^{n/2}\Gamma (\alpha /2)}{\Gamma ((n-\alpha )/2)}.$$ The constant ${\gamma }_n(\alpha )$ is chosen so that $I_{\alpha }{I}_{\beta }=I_{\alpha +\beta }$ whenever $\alpha +\beta <n$, and so that $I_{\alpha }$ is the convolution operator with the fractional kernel $(-\Delta {)}^{-\alpha /2}$ on the Fourier side. The Newtonian potential is $I_2$ in $n\ge 3$ ^{[Riesz 1923]}: explicitly, $I_{2}f(x)=\int f(y)/(|x-y{|}^{n-2}\cdot {\gamma }_n(2))dy$, recovering the Newtonian-potential formula with the correct normalisation.

Theorem 2 (Hardy-Littlewood-Sobolev inequality; Hardy-Littlewood 1928, Sobolev 1938). Let $0<\alpha <n$, and let $1<p<q<\infty$ satisfy $$\frac{1}{q}=\frac{1}{p}-\frac{\alpha }{n}.$$ Then the Riesz potential $I_{\alpha }$ is bounded from $L^p({\mathbb{R}}^n)$ to $L^q({\mathbb{R}}^n)$: $$\Vert I_{\alpha }f{\Vert }_{L^q({\mathbb{R}}^n)}\le C(n,p,\alpha )\Vert f{\Vert }_{L^p({\mathbb{R}}^n)},$$ with $C$ depending only on $n,p,\alpha$ ^{[Hardy-Littlewood 1928] [Sobolev 1938]}. The inequality fails at the endpoints $p=1$ (replaced by a weak-type estimate) and $p=n/\alpha$ (replaced by an $L^{\infty }$-BMO endpoint). Sobolev's contribution was the $\alpha =1$ Sobolev embedding $W^{1,p}({\mathbb{R}}^n)\hookrightarrow L^{p^{\ast }}$ with $1/p^{\ast }=1/p-1/n$; Hardy-Littlewood had proved the one-dimensional inequality $\alpha =1$, $n=1$ in their 1928 paper.

Theorem 3 (Calderón-Zygmund $L^p$ theory; Calderón-Zygmund 1952). Let $f\in L^p({\mathbb{R}}^n)$ for some $1<p<\infty$ and let $u=\Phi \ast f$ be the Newtonian potential (interpreted distributionally). Then for every pair of indices $i,j\in \{1,\dots ,n\}$, the second-order partial derivative ${\partial }_i{\partial }_{j}u$ exists in $L^p({\mathbb{R}}^n)$, and $$\Vert {\partial }_i{\partial }_{j}u{\Vert }_{L^p({\mathbb{R}}^n)}\le C(n,p)\Vert f{\Vert }_{L^p({\mathbb{R}}^n)}\text{for all}1<p<\infty ,$$ where $C(n,p)$ blows up as $p\to 1$ or $p\to \infty$ ^{[Calderón-Zygmund 1952]}. The operator $T_{ij}={\partial }_i{\partial }_j\Phi \ast$ is the prototypical singular integral operator: its convolution kernel ${\partial }_i{\partial }_j\Phi$ is locally integrable only on the boundary of integrability, requiring the principal-value interpretation, and its $L^p$-boundedness is the content of the Calderón-Zygmund theorem.

Theorem 4 (Schauder estimates; Schauder 1934). Let $f\in C_{\mathrm{loc}}^{\alpha }(\Omega )$ for some $0<\alpha <1$ (Hölder-continuous with exponent $\alpha$). If $u\in C^2(\Omega )$ solves $-\Delta u=f$ on $\Omega$, then $u\in C_{\mathrm{loc}}^{2,\alpha }(\Omega )$, and on every compactly contained ${\Omega }^{\prime }⋐\Omega$, $$\Vert u{\Vert }_{C^{2,\alpha }({\Omega }^{\prime })}\le C({\Omega }^{\prime },\Omega ,n,\alpha )\left(\Vert u{\Vert }_{C(\Omega )}+\Vert f{\Vert }_{C^{\alpha }(\Omega )}\right).$$ The Schauder estimates are the Hölder-space companion of the Calderón-Zygmund $L^p$ estimates, and provide the foundational interior regularity for second-order elliptic equations. Their boundary version (Agmon-Douglis-Nirenberg 1959 Comm. Pure Appl. Math. 12) extends the estimate up to the boundary on $C^{2,\alpha }$ domains.

Theorem 5 (Helmholtz decomposition). Every vector field $F\in L^2({\mathbb{R}}^n;{\mathbb{R}}^n)$ has a unique orthogonal decomposition $$F=\nabla u+V,$$ where $u:{\mathbb{R}}^n\to \mathbb{R}$ is a scalar potential (with $\nabla u\in L^2$) and $V:{\mathbb{R}}^n\to {\mathbb{R}}^n$ is divergence-free (i.e. $\mathrm{div}V=0$ distributionally). The scalar potential is given explicitly as $u=-\Phi \ast (\mathrm{div}F)$, recovering the Newtonian potential of the divergence. The Helmholtz decomposition is the foundational tool of vector calculus on ${\mathbb{R}}^n$ and underlies the Hodge decomposition of forms on Riemannian manifolds, the Leray projection in fluid mechanics (Leray 1934 Acta Math. 63), and the gauge-fixing apparatus in electromagnetism and Yang-Mills theory.

Theorem 6 (layer potentials and Fredholm integral equations; Fredholm 1900). Let $\Omega \subset {\mathbb{R}}^n$ be a bounded $C^2$ domain with boundary $\partial \Omega$. For a continuous density $\phi :\partial \Omega \to \mathbb{R}$, define the single-layer potential $S\phi$ and double-layer potential $D\phi$: $$S\phi (x)={\int }_{\partial \Omega }\Phi (x-y)\phi (y)dS(y),D\phi (x)={\int }_{\partial \Omega }{\partial }_{{\nu }_y}\Phi (x-y)\phi (y)dS(y).$$ Both functions are harmonic on $\Omega$ and on ${\mathbb{R}}^n\setminus \overline{\Omega }$. The Dirichlet problem $-\Delta u=0$ on $\Omega$ with $u{|}_{\partial \Omega }=g$ admits a solution of the form $u=D\phi$ in $\Omega$, where $\phi$ satisfies the boundary integral equation $$\frac{1}{2}\phi (x)+{\int }_{\partial \Omega }{\partial }_{{\nu }_y}\Phi (x-y)\phi (y)dS(y)=g(x),x\in \partial \Omega .$$ This is a second-kind Fredholm integral equation, solvable by the Fredholm alternative (Fredholm 1900 Kongl. Vetenskaps-Akademiens Förhandlinger 57) ^{[Fredholm 1900]}: the equation $(I/2+K)\phi =g$ for the compact operator $K$ has a unique solution for every $g$ if and only if the homogeneous equation has only the zero solution, which is a direct application of the maximum principle on $\Omega$.

The Fredholm theory of layer potentials is the historical breakthrough that converted the Dirichlet problem from a delicate existence question into a concrete linear-algebraic problem on the boundary. Its numerical analogue is the boundary element method, where one discretises the boundary integral equation and inverts a finite linear system, with complexity scaling as the boundary discretisation rather than the volume discretisation.

Theorem 7 (Hörmander-Ehrenpreis-Malgrange existence; Ehrenpreis 1954, Malgrange 1955). Every non-zero linear partial differential operator with constant coefficients on ${\mathbb{R}}^n$ admits a fundamental solution. Explicitly, if $P(D)={\sum }_{|\alpha |\le m}{a}_{\alpha }{D}^{\alpha }$ is a non-zero polynomial in the partial-derivative operators $D=({\partial }_1,\dots ,{\partial }_n)$, then there exists a distribution $\Phi \in {\mathcal{D}}^{\prime }({\mathbb{R}}^n)$ with $P(D)\Phi ={\delta }_0$ ^{[Ehrenpreis 1954] [Malgrange 1955]}. Equivalently, for every $f\in C_c^{\infty }({\mathbb{R}}^n)$, the equation $P(D)u=f$ admits a $C^{\infty }$ solution on ${\mathbb{R}}^n$. The result was proved independently and simultaneously by Ehrenpreis (1954 Amer. J. Math. 76) and Malgrange (1955 Ann. Inst. Fourier 6), and is the existence-theoretic capstone of Schwartz's distribution theory (Schwartz 1950 Théorie des distributions).

Theorem 8 (modern applications and computational solvers). The Poisson equation appears in essentially every applied science. Gravitational lensing in general relativity reduces in the weak-field limit to a Poisson equation for the lensing potential, with corrections from the line-of-sight gravitational potential along the photon path. Electrostatic capacitance computation in semiconductor design requires solving Poisson's equation in three dimensions on irregular domains, typically via finite-element discretisation. The pressure equation in computational fluid dynamics ($\Delta p=-\rho \mathrm{div}(u\cdot \nabla u)$ for an incompressible Navier-Stokes flow) is a Poisson equation that must be solved every timestep. Molecular dynamics electrostatics uses Particle-Particle Particle-Mesh (P3M) and Ewald summation, both of which reduce to a Poisson solve on a periodic lattice. Cosmological $N$-body simulations of large-scale structure formation use FFT-based Poisson solvers on the periodic torus, taking advantage of the diagonal action of $-\Delta$ on Fourier modes.

Synthesis. The Poisson equation is the foundational reason that constant-coefficient elliptic theory admits explicit solution operators. The central insight is that the fundamental solution $\Phi$ inverts $-\Delta$ by convolution: every source $f$ with sufficient regularity and decay produces the solution $u=\Phi \ast f$, and this is exactly the integral-operator inverse of the differential operator $-\Delta$. Putting these together with the Riesz-potential calculus of M. Riesz 1923 ^{[Riesz 1923]} and the Hardy-Littlewood-Sobolev inequality of Hardy-Littlewood 1928 ^{[Hardy-Littlewood 1928]} and Sobolev 1938 ^{[Sobolev 1938]}, the convolution-by-$\Phi$ operator extends to a continuous bounded operator $I_2:L^p\to L^q$ on the precise Sobolev-scaling line $1/q=1/p-2/n$. The bridge is between the differential operator $-\Delta$ and its integral inverse $I_2$: the differential side gives the local pointwise PDE, the integral side gives the global representation, and the Sobolev mapping theorem identifies the gain in regularity as exactly two orders.

The pattern generalises through three escalations. First, fractional Riesz potentials $I_{\alpha }=(-\Delta {)}^{-\alpha /2}$ extend the convolution-inverse to non-integer orders and identify the Sobolev spaces $W^{\alpha ,p}$ with the image of $L^p$ under $I_{\alpha }$, the bridge to the modern theory of function spaces. Second, the Calderón-Zygmund framework ^{[Calderón-Zygmund 1952]} generalises the second-derivative-of-Newtonian-potential operator to a class of singular integral operators that remain bounded on $L^p$ for $1<p<\infty$, and identifies the pointwise regularity ${\partial }_i{\partial }_{j}u\in L^p$ when $f\in L^p$, the foundation of modern elliptic $L^p$ theory. Third, the Hörmander-Ehrenpreis-Malgrange existence theorem ^{[Ehrenpreis 1954] [Malgrange 1955]} generalises the Newtonian-potential construction to every non-zero constant-coefficient partial differential operator, identifying $P(D)$ on ${\mathbb{R}}^n$ with an integral operator on Schwartz space whose kernel is the fundamental solution of $P(D)$.

The conceptual closure is the recognition that the Poisson equation is the prototype of every constant-coefficient elliptic boundary-value problem, and its solution apparatus (fundamental solution, Green function, layer potentials, Fredholm integral equations) is the prototype of every linear-elliptic existence and regularity theory. This pattern recurs in the heat equation 02.13.03 via the heat kernel as fundamental solution, in the wave equation via the retarded Green function, in Maxwell's equations via the dyadic Green function, and in general relativity via the Green function of the linearised Einstein operator 13.04.01. The bridge from local PDE to global integral representation, established here for $-\Delta$ in its simplest setting, is the structural fact that makes constant-coefficient PDE theory tractable.

Full proof set Master

Proposition 1 (fundamental solution distributional identity). The fundamental solution $\Phi$ satisfies $-\Delta \Phi ={\delta }_0$ in the sense of distributions on ${\mathbb{R}}^n$. That is, for every test function $\phi \in C_c^{\infty }({\mathbb{R}}^n)$, $$-{\int }_{{\mathbb{R}}^n}\Phi (x)\Delta \phi (x)dx=\phi (0).$$

Proof. Fix $\epsilon >0$ and split the domain into $B_{\epsilon }=B_{\epsilon }(0)$ and its complement: $$-{\int }_{{\mathbb{R}}^n}\Phi \Delta \phi dx=-{\int }_{B_{\epsilon }}\Phi \Delta \phi dx-{\int }_{{\mathbb{R}}^n\setminus B_{\epsilon }}\Phi \Delta \phi dx.$$

For the first integral: $|\Phi |$ is locally integrable (logarithmic in $n=2$, $|x{|}^{2-n}$ in $n\ge 3$), so $|{\int }_{B_{\epsilon }}\Phi \Delta \phi |\le \Vert \Delta \phi {\Vert }_{\infty }{\int }_{B_{\epsilon }}|\Phi |\to 0$ as $\epsilon \to 0$.

For the second integral: on ${\mathbb{R}}^n\setminus B_{\epsilon }$, $\Phi$ is smooth and harmonic. Apply Green's second identity: $${\int }_{{\mathbb{R}}^n\setminus B_{\epsilon }}(\Phi \Delta \phi -\phi \Delta \Phi )dx={\int }_{\partial ({\mathbb{R}}^n\setminus B_{\epsilon })}(\Phi {\partial }_{\nu }\phi -\phi {\partial }_{\nu }\Phi )dS.$$

Since $\Delta \Phi =0$ on ${\mathbb{R}}^n\setminus B_{\epsilon }$ and $\phi$ has compact support (so contributions at infinity vanish), the identity reduces to $${\int }_{{\mathbb{R}}^n\setminus B_{\epsilon }}\Phi \Delta \phi dx=-{\int }_{\partial B_{\epsilon }(0)}(\Phi {\partial }_{\nu }\phi -\phi {\partial }_{\nu }\Phi )dS,$$ where $\nu$ is the unit normal on $\partial B_{\epsilon }$ pointing into $B_{\epsilon }$ (outward from ${\mathbb{R}}^n\setminus B_{\epsilon }$).

The first boundary term has $|\Phi |$ of order ${\epsilon }^{2-n}$ (logarithmic in $n=2$), surface area of order ${\epsilon }^{n-1}$, and $|{\partial }_{\nu }\phi |$ bounded, so $|{\int }_{\partial B_{\epsilon }}\Phi {\partial }_{\nu }\phi dS|=O(\epsilon )$ in $n\ge 3$ (and $O(\epsilon \log \epsilon )$ in $n=2$), tending to zero as $\epsilon \to 0$.

The second boundary term: ${\partial }_{\nu }\Phi (z)=1/(|S^{n-1}|{\epsilon }^{n-1})$ on $\partial B_{\epsilon }(0)$ (computed in Step 5 of the Key Theorem proof). So $${\int }_{\partial B_{\epsilon }(0)}\phi {\partial }_{\nu }\Phi dS=\frac{1}{|S^{n-1}|{\epsilon }^{n-1}}{\int }_{\partial B_{\epsilon }(0)}\phi dS\to \phi (0)$$ by continuity of $\phi$.

Combining and taking $\epsilon \to 0$: $$-{\int }_{{\mathbb{R}}^n}\Phi \Delta \phi dx=\underset{\epsilon \to 0}{\lim }\left(-{\int }_{B_{\epsilon }}\Phi \Delta \phi dx+{\int }_{\partial B_{\epsilon }}\phi {\partial }_{\nu }\Phi dS-{\int }_{\partial B_{\epsilon }}\Phi {\partial }_{\nu }\phi dS\right)=0+\phi (0)-0=\phi (0).$$

$\square$

Proposition 2 (Riesz potential semigroup law). For $0<\alpha ,\beta$ with $\alpha +\beta <n$, the Riesz potentials satisfy $I_{\alpha }{I}_{\beta }=I_{\alpha +\beta }$ on Schwartz space.

Proof. Compute the Fourier transform. The Fourier transform of the Riesz kernel $K_{\alpha }(x)=1/({\gamma }_n(\alpha )|x{|}^{n-\alpha })$ is ${\hat{K}}_{\alpha }(\xi )=|\xi {|}^{-\alpha }$ (the constant ${\gamma }_n(\alpha )$ is chosen precisely to make this so). Therefore $$\widehat{I_{\alpha }f}(\xi )=|\xi {|}^{-\alpha }\hat{f}(\xi ),\widehat{I_{\beta }f}(\xi )=|\xi {|}^{-\beta }\hat{f}(\xi ).$$ Composing: $$\widehat{I_{\alpha }{I}_{\beta }f}(\xi )=|\xi {|}^{-\alpha }\widehat{I_{\beta }f}(\xi )=|\xi {|}^{-\alpha }|\xi {|}^{-\beta }\hat{f}(\xi )=|\xi {|}^{-(\alpha +\beta )}\hat{f}(\xi )=\widehat{I_{\alpha +\beta }f}(\xi ).$$ By Fourier inversion (legitimate on Schwartz space), $I_{\alpha }{I}_{\beta }f=I_{\alpha +\beta }f$ pointwise.

$\square$

Proposition 3 (Hardy-Littlewood-Sobolev inequality, statement and proof outline). For $1<p<q<\infty$ and $0<\alpha <n$ with $1/q=1/p-\alpha /n$, $$\Vert I_{\alpha }f{\Vert }_{L^q({\mathbb{R}}^n)}\le C(n,p,\alpha )\Vert f{\Vert }_{L^p({\mathbb{R}}^n)}.$$

Proof outline. The Marcinkiewicz interpolation theorem reduces the bound to two endpoint estimates: weak-type $(1,n/(n-\alpha ))$ and weak-type $(n/\alpha ,\infty )$. The weak-type $(1,n/(n-\alpha ))$ estimate is the Hardy-Littlewood maximal-function bound for the Riesz kernel.

Concretely: fix $\lambda >0$ and bound $|\{x:|I_{\alpha }f(x)|>\lambda \}|$. Split the kernel into a near-part $K_{\delta }^{<}=K_{\alpha }{1}_{|x|<\delta }$ and a far-part $K_{\delta }^{>}=K_{\alpha }{1}_{|x|\ge \delta }$, and split the convolution accordingly: $I_{\alpha }f=K_{\delta }^{<}\ast f+K_{\delta }^{>}\ast f$. Apply Hölder's inequality to bound the far-part by a constant times $\Vert f{\Vert }_{L^1}^{1-\theta }\Vert f{\Vert }_{L^p}^{\theta }$ for suitable $\theta$, and apply the maximal-function bound to the near-part. Optimising the split parameter $\delta$ gives the weak-type bound.

Marcinkiewicz interpolation then upgrades weak-type to strong-type, recovering the stated $L^p\to L^q$ inequality with the constant $C(n,p,\alpha )$.

The original proofs of Hardy-Littlewood 1928 (for $n=1$) and Sobolev 1938 (for general $n$) used direct rearrangement and integration techniques predating Marcinkiewicz; the modern Calderón-Zygmund decomposition argument is in Stein 1970 Singular Integrals, ch. V.

$\square$

Connections Master

• Laplace equation and harmonic functions 02.13.01. The Poisson equation is the source-driven extension of the Laplace equation. The Newtonian potential construction here is the canonical inversion operator for the Laplacian: every solution of $-\Delta u=f$ on ${\mathbb{R}}^n$ with sufficient decay decomposes as a Newtonian potential of $f$. The Green function on a bounded domain $\Omega$ generalises the construction to allow Dirichlet boundary data, and reduces the solution of $-\Delta u=f$ with $u=g$ on $\partial \Omega$ to the homogeneous-source problem of constructing $G$.

• Lebesgue integral and dominated convergence 02.07.04. Supplies the integration framework on which the Newtonian potential rests. The convolution $\Phi \ast f$ is a Lebesgue integral, and the differentiation-under-the-integral arguments of the Key Theorem proof rely on the dominated convergence theorem to pass derivatives through the integral. The Riesz-potential extension to $L^p$ data, the Calderón-Zygmund $L^p$ theory, and the layer-potential machinery all assume the full Lebesgue-integration apparatus.

• Heat equation 02.13.03. The parabolic analogue of the Poisson equation. The heat kernel ${\Phi }_t(x)=(4\pi t{)}^{-n/2}{e}^{-|x{|}^2/(4t)}$ is the fundamental solution of ${\partial }_t-\Delta$ on ${\mathbb{R}}^n\times (0,\infty )$, and the same Duhamel convolution recipe gives the solution of the inhomogeneous Cauchy problem ${\partial }_{t}u-\Delta u=f$ with $u(\cdot ,0)=u_0$. The steady-state limit ($t\to \infty$ with zero source) recovers harmonic functions, identifying the Poisson equation as the equilibrium of the heat equation.

• Coulomb's law and Gauss's law 10.01.01. The electrostatic version of the Poisson equation. Coulomb's force on a unit charge at position $x$ due to a point charge $q$ at the origin is $-\nabla V(x)$ where $V(x)=q/(4\pi {\epsilon }_0|x|)$, the Newtonian potential of a point source with density $q{\delta }_0/{\epsilon }_0$. Gauss's law $\nabla \cdot E=\rho /{\epsilon }_0$ is the integral form of $-\Delta V=\rho /{\epsilon }_0$, identifying the Poisson equation with the divergence-of-flux statement on charge density. The shell-theorem computation in the Beginner worked example is Newton's gravitational version of the same calculation.

• Einstein field equations 13.04.01. The Poisson equation is the weak-field, slow-motion Newtonian limit of general relativity. The linearised Einstein equations on a Minkowski background reduce, after suitable gauge fixing and dropping time derivatives, to a Poisson equation $-\Delta \Phi =4\pi G\rho$ for the Newtonian gravitational potential, where $\Phi =-h_{00}/2$ in terms of the metric perturbation. Modern post-Newtonian expansions and the parametrised post-Newtonian formalism build on this base, treating relativistic corrections as perturbative additions to the Poisson-equation solution.

Historical & philosophical context Master

Poisson's 1813 note in the Bulletin de la Société Philomathique de Paris ^{[Poisson 1813]} introduced the equation that bears his name, extending Laplace's 1782 equation to allow non-zero sources. Poisson's interest was celestial-mechanics: the gravitational potential inside an attracting body of mass density $\rho$ satisfies $-\Delta \Phi =4\pi G\rho$, whereas outside the body the potential satisfies the source-free Laplace equation. Poisson's note showed how to handle the transition region between exterior and interior potential, opening the modern theory of source-driven potentials.

Gauss's 1840 Allgemeine Lehrsätze ^{[Gauss 1840]} gave the first systematic treatment of what is now called potential theory, including the mean-value theorem, the integration-by-parts identity for the Newtonian potential, and the explicit construction of the gravitational potential of a uniform-density spheroid. Gauss's monograph fixed the modern conventions for the Newtonian-potential normalisation and introduced capacity-theoretic methods (Gauss's theorem on mass attraction) that anticipated the modern Frostman-Choquet theory by a century.

Riesz's 1923 Acta Sci. Math. Szeged paper ^{[Riesz 1923]} introduced fractional integrals as a unifying framework for potential theory. Riesz's $I_{\alpha }$ generalises the Newtonian potential to non-integer orders and identifies the operator $(-\Delta {)}^{-\alpha /2}$ with a convolution against a fractional power of the distance. The Riesz framework was the conceptual breakthrough that converted the Newtonian potential from a single explicit formula into a one-parameter family of operators with a clean composition law $I_{\alpha }{I}_{\beta }=I_{\alpha +\beta }$.

Hardy-Littlewood's 1928 Math. Z. paper ^{[Hardy-Littlewood 1928]} proved the one-dimensional fractional-integral inequality, the prototype of the inequality now bearing the joint Hardy-Littlewood-Sobolev name. Sobolev's 1938 Mat. Sb. paper ^{[Sobolev 1938]} extended the inequality to general $n$ and recognised its connection to the embedding theorems for the spaces $W^{1,p}({\mathbb{R}}^n)$ that Sobolev himself had introduced. The Hardy-Littlewood-Sobolev inequality was the first $L^p$-mapping property of the Newtonian-potential operator and remains the canonical statement of the Sobolev-scaling regularity gain.

Calderón-Zygmund's 1952 Acta Math. paper ^{[Calderón-Zygmund 1952]} established the $L^p$-boundedness of singular integral operators with Calderón-Zygmund kernels, with the second-derivative Newtonian-potential operator ${\partial }_i{\partial }_j\Phi \ast$ as the prototype. The Calderón-Zygmund theorem was the foundational result of modern harmonic analysis and converted the Poisson-equation regularity theory from a $C^k$-Schauder framework to an $L^p$-Sobolev framework. The Calderón-Zygmund machinery (the decomposition lemma, the maximal function, the BMO-John-Nirenberg theory) underlies essentially all subsequent work in elliptic regularity, including the De Giorgi-Nash-Moser theorems for divergence-form equations.

Fredholm's 1900 Kongl. Vetenskaps-Akademiens Förhandlinger paper ^{[Fredholm 1900]} gave the integral-equation reduction of the Dirichlet problem and proved what is now called the Fredholm alternative: a second-kind integral equation $\phi -K\phi =g$ with compact operator $K$ has a unique solution for every $g$ if and only if the homogeneous equation has only the zero solution. Fredholm's theorem converted the Dirichlet problem from a delicate analytic existence question into a concrete linear-algebraic problem on the boundary, and the layer-potential apparatus he introduced remains the canonical foundation of boundary integral methods. The numerical analogue, the boundary element method, is a standard tool in engineering computation.

Schwartz's 1950 Théorie des Distributions ^{[Schwartz 1950]} gave the rigorous foundation for the distributional interpretation of the fundamental solution. Before Schwartz, the identity $-\Delta \Phi ={\delta }_0$ was a formal manipulation in physics texts; after Schwartz, it became a precise statement about a dual-space pairing. The Ehrenpreis-Malgrange existence theorem ^{[Ehrenpreis 1954] [Malgrange 1955]}, proved independently in 1954-55, showed that every non-zero constant-coefficient linear PDE on ${\mathbb{R}}^n$ admits a distributional fundamental solution, generalising the Newtonian-potential construction to the entire class of constant-coefficient operators. Hörmander's subsequent work systematised the existence theory and identified the analytic obstructions to local solvability for variable-coefficient operators.

The Poisson equation now appears as the foundational example in essentially every textbook of partial differential equations and mathematical physics, and the Newtonian-potential apparatus underlies modern computational physics in domains as varied as cosmological $N$-body simulation (FFT-based Poisson solvers on periodic lattices), semiconductor design (finite-element Poisson solvers for electrostatic capacitance), molecular dynamics (Ewald summation and Particle-Particle Particle-Mesh methods), and gravitational lensing (Poisson-equation solvers with general-relativistic corrections). The arc from Poisson's 1813 note to modern teraflop FFT-based solvers is a two-century lineage in which the same equation, the same fundamental solution, and the same convolution recipe have been continuously refined into ever more efficient and ever more general computational tools.

Bibliography Master

@article{Poisson1813,
  author  = {Poisson, Sim\'eon Denis},
  title   = {Remarques sur une \'equation qui se pr\'esente dans la th\'eorie des attractions des sph\'ero\"ides},
  journal = {Bulletin de la Soci\'et\'e Philomathique de Paris},
  volume  = {3},
  year    = {1813},
  pages   = {388--392}
}

@book{Gauss1840,
  author    = {Gauss, Carl Friedrich},
  title     = {Allgemeine Lehrs\"atze in Beziehung auf die im verkehrten Verh\"altnisse des Quadrats der Entfernung wirkenden Anziehungs- und Abstossungskr\"afte},
  publisher = {Weidmannsche Buchhandlung},
  address   = {Leipzig},
  year      = {1840}
}

@article{Fredholm1900,
  author  = {Fredholm, Ivar},
  title   = {Sur une nouvelle m\'ethode pour la r\'esolution du probl\`eme de {D}irichlet},
  journal = {Kongl. Vetenskaps-Akademiens F\"orhandlinger},
  volume  = {57},
  year    = {1900},
  pages   = {39--46}
}

@article{HardyLittlewood1928,
  author  = {Hardy, G. H. and Littlewood, J. E.},
  title   = {Some properties of fractional integrals. {I}},
  journal = {Mathematische Zeitschrift},
  volume  = {27},
  year    = {1928},
  pages   = {565--606}
}

@article{Riesz1923,
  author  = {Riesz, Marcel},
  title   = {Sur les ensembles compacts de fonctions sommables},
  journal = {Acta Scientiarum Mathematicarum (Szeged)},
  volume  = {6},
  year    = {1923},
  pages   = {136--142}
}

@article{Sobolev1938,
  author  = {Sobolev, Sergei L.},
  title   = {On a theorem of functional analysis},
  journal = {Matematicheskii Sbornik (N.S.)},
  volume  = {4},
  year    = {1938},
  pages   = {471--497}
}

@article{CalderonZygmund1952,
  author  = {Calder\'on, Alberto P. and Zygmund, Antoni},
  title   = {On the existence of certain singular integrals},
  journal = {Acta Mathematica},
  volume  = {88},
  year    = {1952},
  pages   = {85--139}
}

@article{Ehrenpreis1954,
  author  = {Ehrenpreis, Leon},
  title   = {Solutions of some problems of division. {I}},
  journal = {American Journal of Mathematics},
  volume  = {76},
  year    = {1954},
  pages   = {883--903}
}

@article{Malgrange1955,
  author  = {Malgrange, Bernard},
  title   = {Existence et approximation des solutions des \'equations aux d\'eriv\'ees partielles et des \'equations de convolution},
  journal = {Annales de l'Institut Fourier},
  volume  = {6},
  year    = {1955},
  pages   = {271--355}
}

@book{Schwartz1950,
  author    = {Schwartz, Laurent},
  title     = {Th\'eorie des Distributions},
  publisher = {Hermann},
  address   = {Paris},
  year      = {1950}
}

@book{Evans2010,
  author    = {Evans, Lawrence C.},
  title     = {Partial Differential Equations},
  edition   = {2},
  publisher = {American Mathematical Society},
  series    = {Graduate Studies in Mathematics},
  volume    = {19},
  year      = {2010}
}

@book{GilbargTrudinger2001,
  author    = {Gilbarg, David and Trudinger, Neil S.},
  title     = {Elliptic Partial Differential Equations of Second Order},
  edition   = {2},
  publisher = {Springer},
  series    = {Grundlehren der mathematischen Wissenschaften},
  volume    = {224},
  year      = {2001}
}

@book{Folland1995,
  author    = {Folland, Gerald B.},
  title     = {Introduction to Partial Differential Equations},
  edition   = {2},
  publisher = {Princeton University Press},
  year      = {1995}
}

@book{Strauss2008,
  author    = {Strauss, Walter A.},
  title     = {Partial Differential Equations: An Introduction},
  edition   = {2},
  publisher = {Wiley},
  year      = {2008}
}

@book{John1982,
  author    = {John, Fritz},
  title     = {Partial Differential Equations},
  edition   = {4},
  publisher = {Springer},
  year      = {1982}
}

@book{Stein1970,
  author    = {Stein, Elias M.},
  title     = {Singular Integrals and Differentiability Properties of Functions},
  publisher = {Princeton University Press},
  year      = {1970}
}