User:Nick Orbeck/PDE
A partial differential equation (or briefly a PDE) is a mathematical equation that involves two or more independent variables, an unknown function (dependent on those variables), and partial derivatives of the unknown function with respect to the independent variables. The order of a partial differential equation is the order of the highest derivative involved. A solution (or a particular solution) to a partial differential equation is a function that solves the equation or, in other words, turns it into an identity when substituted into the equation. A solution is called general if it contains all particular solutions of the equation concerned.
The term exact solution is often used for second- and higher-order nonlinear PDEs to denote a particular solution (see also Preliminary remarks at Second-Order Partial Differential Equations).
Partial differential equations are used to mathematically formulate, and thus aid the solution of, physical and other problems involving functions of several variables, such as the propagation of heat or sound, fluid flow, elasticity, electrostatics, electrodynamics, etc.
First-Order Partial Differential Equations
General Form of First-Order Partial Differential Equation
A first-order partial differential equation with \(n\)independent variables has the general form \[ F\biggl(x_1,x_2,\dots, x_n,w,\frac{\partial w}{\partial x_1}, \frac{\partial w}{\partial x_2},\dots,\frac{\partial w}{\partial x_n}\biggl)=0, \] where \(w=w(x_1,x_2,\dots, x_n)\) is the unknown function and \(F(\dots)\) is a given function.
Quasilinear Equations. Characteristic System. General Solution
General form of first-order quasilinear PDE
A first-order quasilinear partial differential equation with two independent variables has the general form \[\tag{1} f(x,y,w)\frac{\partial w}{\partial x}+g(x,y,w)\frac{\partial w}{\partial y}=h(x,y,w). \]
Such equations are encountered in various applications (continuum mechanics, gas dynamics, hydrodynamics, heat and mass transfer, wave theory, acoustics, multiphase flows, chemical engineering, etc.).
If the functions \(f\), \(g\), and \(h\) are independent of the unknown \(w\), then equation (1) is called linear.
Characteristic system. General solution
The system of ordinary differential equations \[\tag{2} \frac{dx}{f(x,y,w)}=\frac{dy}{g(x,y,w)}=\frac{dw}{h(x,y,w)} \]
is known as the characteristic system of equation (1). Suppose that two independent particular solutions of this system have been found in the form \[\tag{3} u_{1}(x, y, w)=C_{1},\qquad u_{2}(x,y,w)=C_{2}, \] where \(C_1\) and \(C_2\) are arbitrary constants; such particular solutions are known as integrals of system (2). Then the general solution to equation (1) can be written as \[\tag{4} \Phi(u_{1},u_{2})=0, \] where \(\Phi\) is an arbitrary function of two variables. With equation (4) solved for \(u_2\), one often specifies the general solution in the form \(u_2=\Psi(u_1)\), where \(\Psi(u)\) is an arbitrary function of one variable.
Remark. If \(h(x,y,w)\equiv 0\), then \(w=C_2\) can be used as the second integral in (3).
Example. Consider the linear equation \[ \frac{\partial w}{\partial x}+a\frac{\partial w}{\partial y}=b. \] The associated characteristic system of ordinary differential equations \[ \frac{dx}{1}=\frac{dy}{a}=\frac{dw}{b} \] has two integrals \[ y-ax=C_1,\quad \ w-bx=C_2. \] Therefore, the general solution to this PDE can be written as \(w-bx=\Psi(y-ax)\), or \[ w=bx+\Psi(y-ax), \] where \(\Psi(z)\) is an arbitrary function.
Cauchy Problem: Two Formulations. Solving the Cauchy Problem
Generalized Cauchy problem
Generalized Cauchy problem: find a solution \(w=w(x,y)\) to equation (1) satisfying the initial conditions \[\tag{5} x=\varphi_1(\xi ),\quad y=\varphi_2(\xi ),\quad w=\varphi_3(\xi ), \] where \(\xi\) is a parameter \((\alpha\le \xi \le \beta)\) and the \(\varphi_k(\xi)\) are given functions.
Geometric interpretation: find an integral surface of equation (1) passing through the line defined parametrically by equation (5).
Classical Cauchy problem
Classical Cauchy problem: find a solution \(w=w(x,y)\) of equation (1) satisfying the initial condition \[\tag{6} w=\varphi(y)\quad \hbox{at}\quad x=0, \] where \(\varphi(y)\) is a given function.
It is often convenient to represent the classical Cauchy problem as a generalized Cauchy problem by rewriting condition (6) in the parametric form \[ x=0,\quad y=\xi,\quad w=\varphi(\xi). \]
Existence and uniqueness theorem
If the coefficients \(f\), \(g\), and \(h\) of equation (1) and the functions \(\varphi_k\) in (5) are continuously differentiable with respect to each of their arguments and if the inequalities \(f\varphi'_2-g\varphi'_1\not=0\) and \((\varphi'_1)^2+(\varphi'_2)^2\not=0\) hold along the curve (5), then there is a unique solution to the Cauchy problem (in a neighborhood of the curve (5)).
Procedure of solving the Cauchy problem
The procedure for solving the Cauchy problem (1), (5) involves several steps. First, two independent integrals (3) of the characteristic system (2) are determined. Then, to find the constants of integration \(C_1\) and \(C_2\), the initial data (5) must be substituted into the integrals (3) to obtain \[\tag{7} u_{1}\bigl(\varphi_1(\xi),\varphi_2(\xi),\varphi_3(\xi)\bigr)=C_{1},\qquad u_{2}\bigl(\varphi_1(\xi),\varphi_2(\xi),\varphi_3(\xi)\bigr)=C_{2}. \]
Eliminating \(C_1\) and \(C_2\) from (3) and (7) yields \[\tag{8} \begin{array}{rcl} &&u_{1}(x,y,w)=u_{1}\bigl(\varphi_1(\xi),\varphi_2(\xi),\varphi_3(\xi)\bigr),\\ &&u_{2}(x,y,w)=u_{2}\bigl(\varphi_1(\xi),\varphi_2(\xi),\varphi_3(\xi)\bigr). \end{array} \]
Formulas (8) are a parametric form of the solution to the Cauchy problem (1), (5). In some cases, one may succeed in eliminating the parameter \(\xi\) from relations (8), thus obtaining the solution in an explicit form.
In the cases where first integrals (3) of the characteristic system (2) cannot be found using analytical methods, one should employ numerical methods to solve the Cauchy problem (1), (5) (or (1), (6)).
Second-Order Partial Differential Equations
Linear, Semilinear, and Nonlinear Second-Order PDEs
Linear second-order PDEs and their properties. Principle of linear superposition
A second-order linear partial differential equation with two independent variables has the form \[\tag{9} a(x,y)\frac{\partial^2w}{\partial x^2}+ 2b(x,y)\frac{\partial^2w}{\partial x\,\partial y}+ c(x,y)\frac{\partial^2w}{\partial y^2}= \alpha(x,y)\frac{\partial w}{\partial x}+ \beta(x,y)\frac{\partial w}{\partial y}+ \gamma(x, y)w+\delta(x,y). \]
If \(\delta\equiv 0\), equation (9) is a homogeneous linear equation, and if \(\delta\not\equiv 0\), it is a nonhomogeneous linear equation. The functions \(a(x,y)\), \(b(x,y)\), ..., \(\gamma(x,y)\), \(\delta(x,y)\) are called coefficients of equation (9).
Some properties of a homogeneous linear equation (with \(\delta\equiv 0\)):
- A homogeneous linear equation has a particular solution \(w=0\ .\)
- The principle of linear superposition holds; namely, if \(w_1(x,y)\), \(w_2(x,y)\), ..., \(w_n(x,y)\) are particular solutions to homogeneous linear equation, then the function \(A_1w_1(x,y)+A_2w_2(x,y)+\cdots+A_nw_n(x,y),\) where \(A_1\), \(A_2\), ..., \(A_n\) are arbitrary numbers is also an exact solution to that equation.
- Suppose equation (9) has a particular solution \(\tilde w=\tilde w(x,y;\mu)\) that depends on a parameter \(\mu\), and the coefficients of the linear differential equation are independent of \(\mu\) (but can depend on \(x\) and \(y\)). Then, by differentiating \(\tilde w\) with respect to \(\mu\), one obtains other solutions to the equation, \( \frac{\partial\tilde w}{\partial\mu},\quad \frac{\partial^2\tilde w}{\partial\mu^2},\quad \ldots,\quad \frac{\partial^k\tilde w}{\partial\mu^k},\quad \ldots \)
- Let \(\tilde w=\tilde w(x,y;\mu)\) be a particular solution as described in property 3. Multiplying \(\tilde w\) by an arbitrary function \(\varphi(\mu)\) and integrating the resulting expression with respect to \(\mu\) over some interval \([\mu_1,\mu_2]\), one obtains a new function \( \int_{\mu_1}^{\mu_2}\tilde w(x,y;\mu)\varphi(\mu)\,d\mu, \) which is also a solution to the original homogeneous linear equation.
- Suppose the coefficients of the homogeneous linear equation (9) are independent of \(x\ .\) Then: (i) there is a particular solution of the form \(w=e^{\lambda x}u(y)\), where \(\lambda\) is an arbitrary number and \(u(y)\) is determined by a linear second-order ordinary differential equation, and (ii) differentiating any particular solution with respect to \(x\) also results in a particular solution to equation (9).
Properties 2–5 are widely used for constructing solutions to problems governed by linear PDEs.
Examples of particular solutions to linear PDEs can be found in the subsections Heat equation and Laplace equation below.
Semilinear and nonlinear second-order PDEs
A second-order semilinear partial differential equation with two independent variables has the form \[\tag{10} a(x,y)\frac{\partial^2w}{\partial x^2}+ 2b(x,y)\frac{\partial^2w}{\partial x\,\partial y}+ c(x, y)\frac{\partial^2w}{\partial y^2}= F\biggl(x,y,w,\frac{\partial w}{\partial x},\frac{\partial w}{\partial x}\biggr). \]
In the general case, a second-order nonlinear partial differential
equation with two independent variables has the form
\[
F\biggl(x,y,w,\frac{\partial w}{\partial x},\frac{\partial w}{\partial y},
\frac{\partial^2w}{\partial x^2},\frac{\partial^2w}{\partial x\,\partial y},
\frac{\partial^2w}{\partial y^2}\biggr)=0.
\]
The classification and the procedure for reducing linear and semilinear equations of the form (9) and (10) to a canonical form are only determined by the left-hand side of the equations (see below for details).
Some Linear Equations Encountered in Applications
Three basic types of linear partial differential equations are distinguished—parabolic, hyperbolic, and elliptic (for details, see below). The solutions of the equations pertaining to each of the types have their own characteristic qualitative differences.
Heat equation (a parabolic equation)
1. The simplest example of a parabolic equation is the heat equation \[\tag{11} \frac{\partial w}{\partial t}-\frac{\partial^2w}{\partial x^2}=0, \]
where the variables \(t\) and \(x\) play the role of time and a spatial coordinate, respectively. Note that equation (11) contains only one highest derivative term.
Equation (11) is often encountered in the theory of heat and mass transfer. It describes one-dimensional unsteady thermal processes in quiescent media or solids with constant thermal diffusivity. A similar equation is used in studying corresponding one-dimensional unsteady mass-exchange processes with constant diffusivity.
2. The heat equation (11) has infinitely many particular solutions (this fact is common to many PDEs); in particular, it admits solutions \[ \begin{align*} w(x, t)&=A(x^2+2t)+B,\\ w(x, t)&=A\exp(\mu^2t \pm \mu x)+B,\\ w(x, t)&=A\frac 1{\sqrt t}\exp\biggl(-\frac{x^2}{4t}\biggr)+B,\\ w(x, t)&=A\exp(-\mu^2t)\cos(\mu x+B)+C,\\ w(x, t)&=A\exp(-\mu x)\cos(\mu x-2\mu^2t+B)+C, \end{align*} \] where \(A\), \(B\), \(C\), and \(\mu\) are arbitrary constants.
See also Linear heat equations from EqWorld and Heat equation from Wikipedia.
Wave equation (a hyperbolic equation)
1. The simplest example of a hyperbolic equation is the wave equation \[\tag{12} \frac{\partial^2w}{\partial t^2}-\frac{\partial^2w}{\partial x^2}=0, \]
where the variables \(t\) and \(x\) play the role of time and the spatial coordinate, respectively. Note that the highest derivative terms in equation (12) differ in sign.
This equation is also known as the equation of vibration of a string. It is often encountered in elasticity, aerodynamics, acoustics, and electrodynamics.
2. The general solution of the wave equation (12) is \[\tag{13} w=\varphi(x+t)+\psi(x-t), \]
where \(\varphi(x)\) and \(\psi(x)\) are arbitrary twice continuously differentiable functions. This solution has the physical interpretation of two traveling waves of arbitrary shape that propagate to the right and to the left along the \(x\)-axis with a constant speed equal to 1.
See also Wave equation from Wikipedia and Linear hyperbolic equations from EqWorld.
Laplace equation (an elliptic equation)
1. The simplest example of an elliptic equation is the Laplace equation \[\tag{14} \frac{\partial^2w}{\partial x^2}+\frac{\partial^2w}{\partial y^2}=0, \]
where \(x\) and \(y\) play the role of the spatial coordinates. Note that the highest derivative terms in equation (14) have like signs. The Laplace equation is often written briefly as \(\Delta w=0\), where \(\Delta\) is the Laplace operator.
The Laplace equation is often encountered in heat and mass transfer theory, fluid mechanics, elasticity, electrostatics, and other areas of mechanics and physics. For example, in heat and mass transfer theory, this equation describes steady-state temperature distribution in the absence of heat sources and sinks in the domain under study.
A solution to the Laplace equation (14) is called a harmonic function.
2. Note some particular solutions of the Laplace equation (14): \[ \begin{array}{rcl} w(x,y)&=&Ax+By+C,\\ w(x,y)&=&A(x^2-y^2)+Bxy,\\ w(x,y)&=&\displaystyle\frac{Ax+By}{x^2+y^2}+C,\\ w(x,y)&=&(A\sinh\mu x+B\cosh \mu x)(C\cos\mu y+D\sin\mu y),\\ w(x,y)&=&(A\cos \mu x+B\sin \mu x)(C\sinh\mu y+D\cosh\mu y), \end{array} \] where \(A\), \(B\), \(C\), \(D\), and \(\mu\) are arbitrary constants.
A fairly general method for constructing solutions to the Laplace equation (14) involves the following. Let \(f(z)=u(x, y)+iv(x, y)\) be any analytic function of the complex variable \(z=x+iy\) (\(u\) and \(v\) are real functions of the real variables \(x\) and \(y\ ;\) \(i^2=-1\)). Then the real and imaginary parts of \(f\) both satisfy the Laplace equation, \[ \Delta u=0,\qquad \Delta v=0. \] Thus, by specifying analytic functions \(f(z)\) and taking their real and imaginary parts, one obtains various solutions of the Laplace equation (14).
Classification of Second-Order Partial Differential Equations
Types of equations
Any semilinear partial differential equation of the second-order with two independent variables (10) can be reduced, by appropriate manipulations, to a simpler equation that has one of the three highest derivative combinations specified above in examples (11), (12), and (14).
Given a point \((x,y)\), equation (10) is said to be \[ \begin{array}{ll} \mbox{parabolic}& \mbox{if} \ \ b^2-ac=0,\\ \mbox{hyperbolic}&\mbox{if} \ \ b^2-ac>0,\\ \mbox{elliptic}& \mbox{if} \ \ b^2-ac<0 \end{array} \] at this point.
Characteristic equations
In order to reduce equation (10) to a canonical form, one should first write out the characteristic equation \[ a\,(dy)^2-2b\,dx\,dy+c\,(dx)^2=0, \] which with \(a\not\equiv 0\) splits into two equations \[\tag{15} a\,dy-\bigl(b+\sqrt{b^2-ac}\,\bigr)\,dx=0 \]
and \[\tag{16} a\,dy-\bigl(b-\sqrt{b^2-ac}\,\bigr)\,dx=0, \]
and then find their general integrals.
Remark. If \(a\equiv 0\), the simpler equations \[ \begin{array}{rcl} dx&=&0,\\ 2b\,dy-c\,dx&=&0 \end{array} \] should be used instead of (15) and (16). The first equation has the obvious general solution \(x=C\ .\)
Canonical form of parabolic equations (case \(b^2-ac=0\))
In this case, equations (15) and (16) coincide and have a common general integral, \[ u(x, y)=C. \]
By passing from \(x\), \(y\) to new independent variables \(\xi\), \(\eta\) in accordance with the relations \[ \xi=u(x, y),\qquad \eta=\eta(x, y), \] where \(\eta=\eta(x, y)\) is any twice differentiable function that satisfies the condition of nondegeneracy of the Jacobian \(\frac {D(\xi, \eta)}{D(x, y)}\) in the given domain, one reduces equation (10) to the canonical form \[\tag{17} \frac{\partial^2w}{\partial\eta^2}= F_1\biggl(\xi, \eta, w, \frac{\partial w}{\partial \xi}, \frac{\partial w}{\partial \eta}\biggr). \]
As \(\eta\), one can take \(\eta=x\) or \(\eta=y\ .\)
It is apparent that the transformed equation (17) has only one highest-derivative term, just as the heat equation (11).
Two canonical forms of hyperbolic equations (case \(b^2-ac>0\))
1. The general integrals \[ u_1(x, y)=C_1,\qquad u_2(x, y)=C_2 \] of equations (15) and (16) are real and different. These integrals determine two different families of real characteristics.
By passing from \(x\), \(y\) to new independent variables \(\xi\), \(\eta\) in accordance with the relations \[ \xi=u_1(x, y),\qquad \eta=u_2(x, y), \] one reduces equation (10) to \[ \frac{\partial^2w}{\partial\xi\,\partial\eta}= F_2\biggl(\xi, \eta, w, \frac{\partial w}{\partial \xi}, \frac{\partial w}{\partial \eta}\biggr). \] This is the so-called first canonical form of a hyperbolic equation.
2. The transformation \[ \xi=t+z,\qquad \eta=t-z \] brings the above equation to another canonical form, \[ \frac{\partial^2w}{\partial t^2}-\frac{\partial^2w}{\partial z^2}= F_3\biggl(t, z, w, \frac{\partial w}{\partial t}, \frac{\partial w}{\partial z}\biggr), \] where \(F_3=4F_2\ .\) This is the so-called second canonical form of a hyperbolic equation. Apart from notation, the left-hand side of the last equation coincides with that of the wave equation (12).
Canonical form of elliptic equations (case \(b^2-ac<0\))
In this case the general integrals of equations (15) and (16) are complex conjugates; these determine two families of complex characteristics.
Let the general integral of equation (15) have the form \[ u_1(x, y)+iu_2(x, y)=C,\qquad i^2=-1, \] where \(u_1(x, y)\) and \(u_2(x, y)\) are real-valued functions.
By passing from \(x\), \(y\) to new independent variables \(\xi\), \(\eta\) in accordance with the relations \[ \xi=u_1(x, y),\qquad \eta=u_2(x, y), \] one reduces equation (10) to the canonical form \[ \frac{\partial^2w}{\partial\xi^2}+ \frac{\partial^2w}{\partial\eta^2}=F_4\biggl(\xi, \eta, w, \frac{\partial w}{\partial \xi}, \frac{\partial w}{\partial \eta}\biggr). \] Apart from notation, the left-hand side of the last equation coincides with that of the Laplace equation (14).
Basic Problems for PDEs of Mathematical Physics
Most PDEs of mathematical physics govern infinitely many qualitatively similar phenomena or processes. This follows from the fact that differential equations have, as a rule, infinitely many particular solutions. The specific solution that describes the physical phenomenon under study is separated from the set of particular solutions of the given differential equation by means of the initial and boundary conditions.
For simplicity and clarity of illustration, the basic problems of mathematical physics will be presented for the simplest linear equations (11), (12), and (14) only.
Cauchy problem and boundary value problems for parabolic equations
Cauchy problem (\(t\ge 0\), \(-\infty<x<\infty\)). Find a function \(w\) that satisfies heat equation (11) for \(t>0\) and the initial condition \[\tag{18} w=\varphi(x)\quad\hbox{at}\quad t=0. \]
The solution of the Cauchy problem (11), (18) is
\[
w(x, t)=\int^{\infty}_{-\infty}\varphi(\xi)E(x, \xi, t)\,d\xi,
\]
where \(E(x, \xi, t)\) is the fundamental solution of the Cauchy problem,
\[
E(x, \xi, t)=\frac 1{2\sqrt{\pi at}}
\exp\biggl[-\frac{(x-\xi)^2}{4at}\biggr].
\]
In all boundary value problems (or initial-boundary value problems) below, it will be required to find a function \(w\), in a domain \(t\ge 0\ ,\) \(x_1\le x\le x_2\) (\(-\infty<x_1<x_2<\infty\)), that satisfies the heat equation (11) for \(t>0\) and the initial condition (18). In addition, all problems will be supplemented with some boundary conditions as given below.
First boundary value problem. The function \(w(x,t)\) takes prescribed values on the boundary: \[\tag{19} \begin{array}{lll} w=\psi_1(t)& \hbox{at}& x=x_1,\\ w=\psi_2(t)& \hbox{at}& x=x_2. \end{array} \]
In particular, the solution to the first boundary value problem
(11), (18), (19) with
\(\psi_1(t)=\psi_2(t)\equiv 0\), \(x_1=0\), and \(x_2=l\)
is expressed as
\[
w(x,t)=\int^l_0\varphi(\xi)G(x,\xi,t)\,d\xi,
\]
where the Green's function \(G(x,\xi,t)\) is defined by the formulas
\[
\begin{align*}
G(x, \xi, t)&= \frac 2l\sum^{\infty}_{n=1}\sin\biggl(\frac{n\pi x}l\biggr)
\sin\biggl(\frac{n\pi\xi}l\biggr)\exp\biggl(-\frac{an^2\pi^2t}{l^2}\biggr)\\
&=\frac 1{2\sqrt{\pi at}}
\sum^{\infty}_{n=-\infty}\biggl\{\exp\biggl[-\frac{(x-\xi+2nl)^2}{4at}\biggr]-
\exp\biggl[-\frac{(x+\xi+2nl)^2}{4at}\biggr]\biggr\}.
\end{align*}
\]
The first series converges rapidly at large \(t\) and the second series at
small \(t\ .\)
Second boundary value problem. The derivatives of the function \(w(x,t)\) are prescribed on the boundary: \[\tag{20} \begin{aligned} \frac{\partial w}{\partial x}=\psi_1(t)&\quad \hbox{at}\quad x=x_1,\\ \frac{\partial w}{\partial x}=\psi_2(t)&\quad \hbox{at}\quad x=x_2. \end{aligned} \]
Third boundary value problem. A linear relationship between the
unknown function and its derivatives are prescribed on the boundary:
\[\tag{21}
\begin{aligned}
\frac{\partial w}{\partial x}-k_1w=\psi_1(t)&\quad \hbox{at}\quad x=x_1,\\
\frac{\partial w}{\partial x}+k_2w=\psi_2(t)&\quad \hbox{at}\quad x=x_2.
\end{aligned}
\]
Mixed boundary value problems. Conditions of different type, listed
above, are set on the boundary of the domain in question, for example,
\[\tag{22}
\begin{array}{rll}
x=\psi_1(t)& \hbox{at}& x=x_1,\\
\displaystyle\frac{\partial w}{\partial x}=\psi_2(t)& \hbox{at}& x=x_2.
\end{array}
\]
The boundary conditions (19)–(22) are called homogeneous if
\(\psi_1(t)=\psi_2(t)\equiv 0\ .\)
Solutions to the above initial-boundary value problems for the heat equation can be obtained by separation of variables (Fourier method) in the form of infinite series or by the method of integral transforms using the Laplace transform.
For other linear heat equations, their exact solutions, and solutions to associated Cauchy problems and boundary value problems, see Linear heat equations at EqWorld.
Cauchy problem and boundary value problems for hyperbolic equations
Cauchy problem (\(t\ge 0\), \(-\infty<x<\infty\)). Find a function \(w\) that satisfies the wave equation (12) for \(t>0\) and two initial conditions \[\tag{23} \begin{array}{rll} w=\varphi_0(x)& \hbox{at}& t=0,\\ \displaystyle\frac{\partial w}{\partial t}=\varphi_1(x)& \hbox{at}& t=0. \end{array} \]
The solution of the Cauchy problem (12), (23) is given by d'Alembert's formula:
\[
w(x, t)=\frac 12[\varphi_0(x+at)+\varphi_0(x-at)]+\frac
1{2a}\int^{x+at}_{x-at}\varphi_1(\xi)\,d\xi.
\]
Boundary value problems. In all boundary value problems, it is required to find a function \(w\), in a domain \(t\ge 0\), \(x_1\le x\le x_2\) (\(-\infty<x_1<x_2<\infty\)), that satisfies the wave equation (12) for \(t>0\) and the initial conditions (23). In addition, appropriate boundary conditions, (19), (20), (21), or (22), are imposed.
Solutions to these boundary value problems for the wave equation can be obtained by separation of variables (Fourier method) in the form of infinite series. In particular, the solution to the first boundary value problem (12), (19), (23) with homogeneous boundary conditions, \(\psi_1(t)=\psi_2(t)\equiv 0\) at \(x_1=0\) and \(x_2=l\), is expressed as \[\tag{24} w(x,t)=\frac{\partial}{\partial t}\int^l_0\varphi_0(\xi)G(x,\xi,t)\,d\xi +\int^l_0\varphi_1(\xi)G(x,\xi,t)\,d\xi, \]
where \[ G(x, \xi, t) =\frac 2{a\pi}\sum^{\infty}_{n=1}\frac 1n\sin\Bigl(\frac{n\pi x}l\Bigr) \sin\Bigl(\frac{n\pi\xi}l\Bigr)\sin\Bigl(\frac{n\pi at}l\Bigr). \]
Goursat problem. On the characteristics of the wave equation (12), values of the unknown function \(w\) are prescribed: \[\tag{25} \begin{array}{llll} w=\varphi(x)& \hbox{for}& x-t=0& (0\le x\le a),\\ w=\psi(x)& \hbox{for}& x+t=0& (0\le x\le b), \end{array} \]
with the consistency condition \(\varphi(0)=\psi(0)\) implied to hold.
Substituting the values set on the characteristics (25) into the general solution of the wave equation (13), one arrives at a system of linear algebraic equations for \(\varphi(x)\) and \(\psi(x)\ .\) As a result, the solution to the Goursat problem (12), (25) is obtained in the form \[ w(x,t)=\varphi\biggl(\frac{x+t}2\biggr)+\psi\biggl(\frac{x-t}2\biggr)-\varphi(0). \] The solution propagation domain is the parallelogram bounded by the four lines \[ x-t=0,\quad x+t=0,\quad x-t=2b,\quad x+t=2a. \]
For other linear wave equations, their exact solutions, and solutions to associated Cauchy problems and boundary value problems, see Linear hyperbolic equations at EqWorld.
Boundary value problems for elliptic equations
Setting boundary conditions for the first, second, and third boundary value problems for the Laplace equation (14) means prescribing values of the unknown function, its first derivative, and a linear combination of the unknown function and its derivative, respectively.
For example, the first boundary value problem in a rectangular domain \(0\le x\le a\), \(0\le y\le b\) is characterized by the boundary conditions \[\tag{26} \begin{array}{llllll} w=\varphi_1(y)& \hbox{at}& x=0,\quad& w=\varphi_2(y)& \hbox{at}& x=a,\\ w=\varphi_3(x)& \hbox{at}& y=0,\quad& w=\varphi_4(x)& \hbox{at}& y=b. \end{array} \]
The solution to problem (14), (26) with \(\varphi_3(x)=\varphi_4(x)\equiv 0\) is given by
\[
w(x, y)=\sum^{\infty}_{n=1}
A_n\sinh\biggl[\frac{n\pi}b(a-x)\biggr]\sin\biggl(\frac{n\pi}by\biggr)
+\sum^{\infty}_{n=1}
B_n\sinh\biggl(\frac{n\pi}bx\biggr)\sin\biggl(\frac{n\pi}by\biggr),
\]
where the coefficients \(A_n\) and \(B_n\) are expressed as
\[
A_n=\frac{2}{\lambda_n}\int^b_0\varphi_1(\xi)\sin\biggl(\frac{n\pi\xi}b\biggr)d\xi,\quad
B_n=\frac{2}{\lambda_n}\int^b_0\varphi_2(\xi)\sin\biggl(\frac{n\pi\xi}b\biggr)d\xi,\quad
\lambda_n=b\sinh\biggl(\frac{n\pi a}b\biggr).
\]
Remark. For elliptic equations, the first boundary value problem is often called the Dirichlet problem, and the second boundary value problem is called the Neumann problem.
For other linear elliptic equations, their exact solutions, and solutions to associated boundary value problems, see Linear elliptic equations at EqWorld.
Some Nonlinear Equations Encountered in Applications
Nonlinear heat equation:
\[\tag{27} \frac{\partial w}{\partial t}=\frac{\partial}{\partial x} \biggl[f(w)\frac{\partial w}{\partial x}\biggr]. \]
This equation describes one-dimensional unsteady thermal processes in quiescent media or solids in the case where the thermal diffusivity is temperature dependent, \(f(w)>0\ .\) In the special case \(f(w)\equiv 1\), the nonlinear equation (27) becomes the linear heat equation (11).
In general, the nonlinear heat equation (27) admits exact solutions of the form \[ \begin{array}{ll} w=W(kx-\lambda t)& (\hbox{traveling-wave solution}),\\ w=U(x/\!\sqrt t\,)& (\hbox{self-similar solution}), \end{array} \] where \(W=W(z)\) and \(U=U(r)\) are determined by ordinary differential equations, and \(k\) and \(\lambda\) are arbitrary constants.
Kolmogorov–Petrovskii–Piskunov equation:
\[\tag{28} \frac{\partial w}{\partial t}=a\frac{\partial^2w}{\partial x^2}+f(w),\qquad a>0. \]
Equations of this form are often encountered in various problems of mass and heat transfer (with \(f\) being the rate of a volume chemical reaction), combustion theory, biology, and ecology.
In the special case of \(f(w)\equiv 0\) and \(a=1\), the nonlinear equation (28) becomes the linear heat equation (11).
Remark. Equation (28) is also called a heat equation with a nonlinear source.
Burgers equation:
\[\tag{29} \frac{\partial w}{\partial t}+w\frac{\partial w}{\partial x}=\frac{\partial^2w}{\partial x^2}. \]
This equation is used for describing wave processes in gas dynamics, hydrodynamics, and acoustics.
1. Exact solutions to the Burgers equation can be obtained using the following formula (Hopf–Cole transformation): \[ w(x,t)=-\frac 2u\frac{\partial u}{\partial x}, \] where \(u=u(x,t)\) is a solution to the linear heat equation \(u_t=u_{xx}\) (see above for details).
2. The solution to the Cauchy problem for the Burgers equation with the initial condition \[ w=f(x)\quad {\rm at}\quad t=0 \qquad (-\infty<x<\infty) \] has the form \[ w(x,t)=-2\frac {\partial}{\partial x}\ln F(x,t), \] where \[ F(x, t)=\frac 1{\sqrt{4\pi t}}\int^{\infty}_{-\infty} \exp\biggl[-\frac{(x-\xi)^2}{4t}+\frac12\int^{\xi}_0f(\xi')\,d\xi'\biggr]d\xi. \]
Nonlinear wave equation:
\[\tag{30} \frac{\partial^2w}{\partial t^2}=\frac{\partial}{\partial x} \biggl[f(w)\frac{\partial w}{\partial x}\biggr]. \]
This equation is encountered in wave and gas dynamics, \(f(w)>0\ .\) In the special case \(f(w)\equiv 1\), the nonlinear equation (30) becomes the linear wave equation (12).
Equation (30) admits exact solutions in implicit form: \[ \begin{array}{rcl} x+t\sqrt{f(w)}&=&\varphi(w),\\ x-t\sqrt{f(w)}&=&\psi(w), \end{array} \] where \(\varphi(w)\) and \(\psi(w)\) are arbitrary functions.
Equation (30) can be reduced to a linear equation (see Polyanin and Zaitsev, 2004).
Nonlinear Klein–Gordon equation:
\[\tag{31} \frac{\partial^2w}{\partial t^2}=a\frac{\partial^2w}{\partial x^2}+f(w), \qquad a>0. \]
Equations of this form arise in differential geometry and various areas of physics (superconductivity, dislocations in crystals, waves in ferromagnetic materials, laser pulses in two-phase media, and others). For \(f(w)\equiv 0\) and \(a=1\), equation (31) coincides with the linear wave equation (12).
1. In general, the nonlinear Klein–Gordon equation (31) admits exact solutions of the form \[ \begin{array}{ll} w=W(z),& z=kx-\lambda t,\\ w=U(\xi),& \xi=(\sqrt a\,t+C_1)^2-(x+C_2)^2, \end{array} \] where \(W=W(z)\) and \(U=U(\xi)\) are determined by ordinary differential equations, while \(k\), \(\lambda\), \(C_1\), and \(C_2\) are arbitrary constants.
2. In the special case \[ f(w)=be^{\beta w}, \] the general solution of equation (31) is expressed as \[ w(x,t)=\frac 1{\beta}\bigl[\varphi(z)+\psi(y)\bigr]- \frac 2{\beta}\ln\biggl|k\int \exp\bigl[\varphi(z)\bigr]\,dz -\frac{b\beta}{8ak}\int\exp\bigl[\psi(y)\bigr]\,dy\biggr|, \] \[ z=x-\sqrt a\,t,\qquad y=x+\sqrt a\,t, \] where \(\varphi=\varphi(z)\) and \(\psi=\psi(y)\) are arbitrary functions and \(k\) is an arbitrary constant.
Remark. In the special cases \(f(w)=b\sin(\beta w)\) and \(f(w)=b\sinh(\beta w)\), equation (31) is called the sine-Gordon equation and the sinh-Gordon equation, respectively.
Nonlinear Laplace equation:
\[\tag{32} \frac{\partial^2w}{\partial x^2}+\frac{\partial^2w}{\partial y^2}=f(w). \]
This equation is also called a stationary heat equation with a nonlinear source.
1. In general, the nonlinear heat equation (32) admits exact solutions of the form \[ \begin{array}{ll} w=W(z),& z=k_1x+k_2y,\\ w=U(r),& r=\sqrt{(x+C_1)^2+(y+C_2)^2}, \end{array} \] where \(W=W(z)\) and \(U=U(r)\) are determined by ordinary differential equations, while \(k_1\), \(k_2\), \(C_1\), and \(C_2\) are arbitrary constants.
2. In the special case \[ f(w)=ae^{\beta w}, \] the general solution of equation (32) is expressed as \[ w(x,y)=-\frac 2\beta\ln\frac{\bigl|1-2a\beta\Phi(z)\overline{\Phi(z)}\,\bigr|}{4|\Phi'_z(z)|}, \] where \(\Phi=\Phi(z)\) is an arbitrary analytic function of the complex variable \(z=x+iy\) with nonzero derivative, and the bar over a symbol denotes the complex conjugate.
Monge–Ampere equation:
\[ \biggl(\frac{\partial^2w}{\partial x\,\partial y}\biggr)^{\!2}- \frac{\partial^2w}{\partial x^2} \frac{\partial^2w}{\partial y^2}=f(x,y). \] The equation is encountered in differential geometry, gas dynamics, and meteorology.
Below are solutions to the homogeneous Monge–Ampere equation for the special case \(f(x,y)\equiv 0\ .\)
1. Exact solutions involving one arbitrary function: \[ w(x,y)=\varphi(C_1x+C_2 y)+C_3x+C_4y+C_5, \] \[ w(x,y)=(C_1x+C_2y)\,\varphi\biggl(\frac yx\biggr)+C_3x+C_4y+C_5, \] \[ w(x,y)=(C_1x+C_2y+C_3)\,\varphi\biggl(\frac{C_4x+C_5y+C_6}{C_1x+C_2y+C_3}\biggr) +C_7x+C_8y+C_9, \] where \(C_1\), ..., \(C_{9}\) are arbitrary constants and \(\varphi=\varphi(z)\) is an arbitrary function.
2. General solution in parametric form: \[ w=t x+\varphi(t)y+\psi(t), \] \[ x+\varphi'(t)y+\psi'(t)=0, \] where \(t\) is the parameter, and \(\varphi=\varphi(t)\) and \(\psi=\psi(t)\) are arbitrary functions.
Simplest Types of Exact Solutions of Nonlinear PDEs
Preliminary remarks
The following classes of solutions are usually regarded as exact solutions to nonlinear partial differential equations of mathematical physics:
- Solutions expressible in terms of elementary functions.
- Solutions expressed by quadrature.
- Solutions described by ordinary differential equations (or systems of ordinary differential equations).
- Solutions expressible in terms of solutions to linear partial differential equations (and/or solutions to linear integral equations).
The simplest types of exact solutions to nonlinear PDEs are traveling-wave solutions and self-similar solutions. They often occur in various applications.
In what follows, it is assumed that the unknown \(w\) depends on two variables, \(x\) and \(t\), where \(t\) plays the role of time and \(x\) is a spatial coordinate.
Traveling-wave solutions
Traveling-wave solutions, by definition, are of the form \[\tag{33} w(x,t)=W(z),\quad \ z=kx-\lambda t, \]
where \(\lambda/k\) plays the role of the wave propagation velocity (the value \(\lambda =0\) corresponds to a stationary solution, and the value \(k=0\) corresponds to a space-homogeneous solution). Traveling-wave solutions are characterized by the fact that the profiles of these solutions at different time instants are obtained from one another by appropriate shifts (translations) along the \(x\)-axis. Consequently, a Cartesian coordinate system moving with a constant speed can be introduced in which the profile of the desired quantity is stationary. For \(k>0\) and \(\lambda>0\), the wave (33) travels along the \(x\)-axis to the right (in the direction of increasing \(x\)).
Traveling-wave solutions occur for equations that do not explicitly involve independent variables, \[\tag{34} F\biggl(w, \frac{\partial w}{\partial x}, \frac{\partial w}{\partial t}, \frac{\partial^2w}{\partial x^2}, \frac{\partial^2w}{\partial x\,\partial t}, \frac{\partial^2w}{\partial t^2},\ldots\biggr)=0. \]
Substituting (33) into (34), one obtains an autonomous ordinary differential equation for the function \(W(z)\ :\) \[ F(W,kW',-\lambda W',k^2W'',-k\lambda W'',\lambda ^2W'',\ldots)=0, \] where \(k\) and \(\lambda \) are arbitrary constants, and the prime denotes a derivative with respect to \(z\ .\)
Remark. The term traveling-wave solution is also used in the cases where the variable \(t\) plays the role of a spatial coordinate, \(t=y\ .\)
All nonlinear equations considered above, (27)–(32) and (33) with \(f(x,y)=0\), admit traveling-wave solutions.
Self-similar solutions
By definition, a self-similar solution is a solution of the form \[\tag{35} w(x,t)=t^{\alpha}U(\zeta),\quad \ \zeta=xt^\beta. \]
The profiles of these solutions at different time instants are obtained from one another by a similarity transformation (like scaling).
Self-similar solutions exist if the scaling of the independent and dependent variables, \[\tag{36} t=C\bar t,\quad x=C^k\bar x,\quad w=C^m\bar w,\qquad \mbox{where}\ C\not=0\ \mbox{is an arbitrary constant}, \]
for some \(k\) and \(m\) such that \(|k|+|m|\not=0\), is equivalent to the identical transformation.
It can be shown that the parameters in solution (35) and transformation (36) are linked by the simple relations \[\tag{37} \alpha=m, \quad \ \beta=-k. \]
In practice, the above existence criterion is checked and if a pair of \(k\)
and \(m\) in (36) has been found, then a self-similar solution is defined by formulas (35) with parameters (37).
Example. Consider the heat equation with a nonlinear power-law source term \[\tag{38} \frac{\partial w}{\partial t}=a\frac{\partial^2w}{\partial x^2}+bw^n. \]
The scaling transformation (36) converts equation (38) into \[\tag{39} C^{m-1}\frac{\partial \bar w}{\partial\bar t}= aC^{m-2k}\frac{\partial^2\bar w}{\partial \bar x^2}+bC^{mn}\bar w^n. \]
In order that equation (39) coincides with (38), one must require that the powers of \(C\) are the same, which yields the following system of linear algebraic equations for the constants \(k\) and \(m\ :\) \[ m-1=m-2k=mn. \] This system admits a unique solution\[\,k=\frac 12\ ,\] \(m=\frac 1{1-n}\ .\) Using this solution together with relations (35) and (37), one obtains self-similar variables in the form \[ w=t^{1/(1-n)}U(\zeta),\quad \ \zeta=xt^{-1/2}. \] Inserting these into (38), one arrives at the following ordinary differential equation for \(U(\zeta)\ :\) \[ aU''_{\zeta\zeta}+\frac12\zeta U'_\zeta+\frac 1{n-1}U+bU^n=0. \]
Cauchy Problem and Boundary Value Problems for Nonlinear Equations
The Cauchy problem and boundary value problems for nonlinear equations are stated in exactly the same way as for linear equations (see Basic Problems for PDEs of Mathematical Physics).
Examples. The Cauchy problem for a nonlinear heat equation is stated as follows: find a solution to equation (27) subject to the initial condition (18).
The first boundary value problem for a nonlinear wave equation as follows: find a solution to equation (32) subject to the initial conditions (18) and the boundary conditions (19).
Problems for nonlinear PDEs are normally solved using numerical methods.
Higher-Order Partial Differential Equations
Apart from second-order PDEs, higher-order equations also quite often arise in applications. Below are only a few important examples of such equations with some of their solutions.
Higher-Order Linear Partial Differential Equations
Equation of transverse vibration of elastic rod:
\[ \frac{\partial^2w}{\partial t^2}+a^2\frac{\partial^4w}{\partial x^4}=0. \]
The equation has the following particular solutions: \[ \begin{array}{l} w(x,t)=\bigl[A\sin(\lambda x)+B\cos(\lambda x)+C\sinh(\lambda x)+D\cos(\lambda x)\bigr]\sin(\lambda^2at),\\ w(x,t)=\bigl[A_1\sin(\lambda x)+B_1\cos(\lambda x)+C_1\sinh(\lambda x)+ D_1\cos(\lambda x)\bigr]\cos(\lambda^2at), \end{array} \] where \(A\), \(B\), \(C\), \(D\ ,\) \(A_1\), \(B_1\), \(C_1\), \(D_1\), and \(\lambda\) are arbitrary constants.
For solutions to associated Cauchy problems and boundary value problems, see Equation of transverse vibration of elastic rods at EqWorld.
Biharmonic equation:
\[\tag{40} \Delta\Delta w=0, \]
where \(\Delta\Delta\) is the biharmonic operator, \[ \Delta\Delta \equiv\Delta^{\!2}= \frac{\partial^4}{\partial x^4}+ 2\frac{\partial^4}{\partial x^2\,\partial y^2}+ \frac{\partial^4}{\partial y^4}. \]
The biharmonic equation (40) is encountered in plane problems of elasticity (\(w\)is the Airy stress function). It is also used to describe slow flows of viscous incompressible fluids (\(w\) is the stream function).
Various representations of the general solution to equation (40) in terms of harmonic functions include \[ \!\!\!\begin{array}{l} w(x,y)=xu_1(x,y)+u_2(x,y),\\ w(x,y)=yu_1(x,y)+u_2(x,y),\\ w(x,y)=(x^2+y^2)u_1(x,y)+u_2(x,y), \end{array} \] where \(u_1\) and \(u_2\) are arbitrary functions satisfying the Laplace equation \(\Delta u_k=0\,\) (\(k=1,\,2\)).
Complex form of representation of the general solution: \[ w(x,y)=\mbox{Re}\bigl[\overline z f(z)+g(z)\bigr], \] where \(f(z)\) and \(g(z)\) are arbitrary analytic functions of the complex variable \(z=x+iy\ ;\) \(\overline z=x-iy\), \(i^2=-1\ .\) The symbol \(\mbox{Re}[A]\) stands for the real part of a complex quantity \(A\ .\)
For solutions to associated boundary value problems, see Biharmonic equation at EqWorld.
Higher-Order Nonlinear Partial Differential Equations
Korteweg–de Vries equation:
\[ \frac{\partial w}{\partial t}+\frac{\partial^3w}{\partial x^3} -6w\frac{\partial w}{\partial x}=0. \] It is used in many sections of nonlinear mechanics and theoretical physics for describing one-dimensional nonlinear dispersive nondissipative waves. In particular, the mathematical modeling of moderate-amplitude shallow-water surface waves is based on this equation. For exact solutions to this equation, see Korteweg–de Vries equation at EqWorld.
Equation of a steady laminar boundary layer on a flat plate:
\[ \frac{\partial w}{\partial y}\frac{\partial^2w}{\partial x\,\partial y}- \frac{\partial w}{\partial x}\frac{\partial^2w}{\partial y^2}=a \frac{\partial^3w}{\partial y^3}. \] where \(w\) is the stream function. For exact solutions, see Boundary layer equations at EqWorld.
Boussinesq equation:
\[ \frac{\partial^2w}{\partial t^2}+\frac{\partial}{\partial x} \biggl(w\frac{\partial w}{\partial x}\biggr)+\frac{\partial^4w}{\partial x^4}=0. \] This equation arises in several physical applications: propagation of long waves in shallow water, one-dimensional nonlinear lattice-waves, vibrations in a nonlinear string, and ion sound waves in a plasma. For exact solutions, see Boussinesq equation at EqWorld.
Equation of motion of a viscous fluid:
\[ \frac{\partial w}{\partial y}\frac{\partial}{\partial x}(\Delta w)- \frac{\partial w}{\partial x}\frac{\partial}{\partial y}(\Delta w)=a\,\Delta\Delta w,\qquad \Delta w=\frac{\partial^2w}{\partial x^2}+\frac{\partial^2w}{\partial y^2}. \] This is a two-dimensional stationary equation of motion of a viscous incompressible fluid—it is obtained from the Navier–Stokes equations by the introduction of the stream function \(w\ .\) For exact solutions to this equation, see Navier–Stokes equations at EqWorld.
Approximate and Numerical Methods
The preceding discussion pertains to the exact or analytical solution of PDEs. For example, in the case of a heat equation or a wave equation, an exact solution would be a function \(w=f(x,t)\) which, when substituted into the respective equation would satisfy it identically along with all of the associated initial and boundary conditions.
Although analytical solutions are exact, they also may not be available, simply because we do not know how to derive such solutions. This could be because the PDE system has too many PDEs, or they are too complicated, e.g., nonlinear, or both, to be amenable to analytical solution. In this case, we may have to resort to an approximate solution. That is, we seek an analytical or numerical approximation to the exact solution.
Perturbation methods are an important subset of approximate analytical methods. They may be applied if the problem involves small (or large) parameters, which are used for constructing solutions in the form of asymptotic expansions. For books on perturbation methods, see Google Book Search. These and other methods for PDEs are also outlined in Zwillinger (1997).
Unlike exact and approximate analytical methods, methods to compute numerical PDE solutions are in principle not limited by the number or complexity of the PDEs. This generality combined with the availability of high performance computers makes the calculation of numerical solutions feasible for a broad spectrum of PDEs (such as the Navier–Stokes equations) that are beyond analysis by analytical methods. The development and implementation (as computer codes) of numerical methods or algorithms for PDE systems is a very active area of research. Here we indicate in the external links just two readily available links to Scholarpedia.
We now consider the numerical solution of a parabolic PDE, a hyperbolic PDE, and an elliptic PDE.
Parabolic PDE
Analytical solutions to a parabolic PDE (heat equation) are given here. But we will proceed with a numerical solution and use one of these analytical solutions to evaluate the numerical solution.
We can consider the numerical solution to the heat equation \[\tag{41} \frac{\partial w}{\partial t}-\frac{\partial^2w}{\partial x^2}=0 \]
as a two-step process:
- Numerical approximation of the derivative \(\dfrac{\partial ^2 w}{\partial x^2}\ .\) At this point, we will have a semi-discretization of Eq. (41).
- Numerical approximation of the derivative \(\dfrac{\partial w}{\partial t}\ .\) At this point, we will have a full discretization of Eq. (41).
In order to implement these two steps, we require a grid in \(x\) and a grid in \(t\ .\) For the grid in \(x\), we denote a position along the grid with the index \(i\ .\) Then we can consider the Taylor series expansion of the numerical solution at grid point \(i\) \[\tag{42} w_{i+1}=w_{i} + \dfrac{dw_i}{dx}(x_{i+1}-x_i) + \dfrac{d^2w_i}{dx^2}\dfrac{(x_{i+1}-x_{i})^2}{2!} + \dfrac{d^3w_i}{dx^3}\dfrac{(x_{i+1}-x_{i})^3}{3!} + \cdots \]
and \[\tag{43} w_{i-1}=w_{i} + \dfrac{dw_i}{dx}(x_{i-1}-x_i) + \dfrac{d^2w_i}{dx^2}\dfrac{(x_{i-1}-x_{i})^2}{2!} + \dfrac{d^3w_i}{dx^3}\dfrac{(x_{i-1}-x_{i})^3}{3!} + \cdots \]
If we consider a uniform grid (a grid with uniform spacing \(\Delta x = x_{i+1}-x_i = x_{i} - x_{i-1})\), addition of Eqs. (42) and (43) gives (note the cancellation of the first and third derivative terms since \(x_{i-1}-x_{i} = -\Delta x\)) \[\tag{44} w_{i+1}+w_{i-1}=2w_{i} + \dfrac{d^2w_i}{dx^2}\Delta x^2 + O(\Delta x^4) \]
where \(O(\Delta x^4)\) denotes a term proportional to \(\Delta x^4\) or of order \(\Delta x^4\ ;\) this term can be considered a truncation error resulting from truncating the Taylor series of Eqs. (42) and (43) beyond the \(\Delta x^2\) term. Then Eq. (44) gives for the second derivative \[\tag{45} \dfrac{d^2w_i}{dx^2} \approx \dfrac{w_{i+1}-2w_{i}+w_{i-1}}{\Delta x^2} + O(\Delta x^2) \]
Equation (45) is a second order (because of the principal error or truncation error \(O(\Delta x^2)\)) finite difference approximation of \(d^2w_i/dx^2\ .\)
If Eq. (45) is substituted in Eq. (41) (to replace the derivative \(\frac{\partial ^2w} {\partial x^2}\)), a system of ODEs results \[\tag{46} \dfrac{dw_i}{dt}=D\dfrac{w_{i+1}-2w_{i}+w_{i-1}}{\Delta x^2}+O(\Delta x^2),\quad i=1,2,\dots,N \]
(we have added a multiplying constant \(D\) to the right-hand side of Eq. (41), generally termed a thermal diffusivity if \(w\) in Eq. (41) is temperature and a mass diffusivity if \(w\) is concentration; \(D\) has the MKS units m2/s as expected from a consideration of Eq. (41) with \(x\) in metres and \(t\) in seconds).
Note that the independent variable \(x\) does not appear explicitly in Eqs. (46) and that the only independent variable is \(t\) (so that they are ODEs). \(N\) is the number of points in the \(x\) grid (\(x\) is termed a boundary value variable since the terminal grid points at \(i=1\) and \(i=N\) typically refer to the boundaries of a physical system). Thus Eq. (41) is partly discretized (in \(x\)) and therefore Eqs. (46) are referred to as a semi-discretization.
To compute a solution to Eq. (41), we could apply an established initial-value integrator in \(t\ .\) This is the essence of the method of lines (MOL). Alternatively, we could now discretize Eqs. (46). For example, if we apply Eq. (42) on a grid in \(t\) with an index \(k\) \[\tag{47} w_{i}^{k+1}=w_i^k + \dfrac{dw_i^k}{dt}(t^{k+1}-t^k) + \dfrac{d^2w_i^k}{dt^2} \dfrac{(t^{k+1}-t^{k})^2}{2!} + \cdots,\quad i=1,2,\dots,N,\quad k=1,2,\dots \]
If the grid in \(t\) has a uniform spacing \(h = t^{k+1}-t^k\) and if truncation after the first derivative term is applied, \[\tag{48} w_{i}^{k+1}=w_i^k + \dfrac{dw_i^k}{dt}h + O(h^2) \]
Equation (48), the classical Euler's method, can be used to step along the solution of Eq. (41) from point \(k\) to \(k+1\) (at a grid point \(i\) in \(x\)).
Application of Eq. (48) to Eq. (46) gives the fully discretized approximation of Eq. (41) \[\tag{49} w_i^{k+1}=w_i^k+hD\dfrac{w_{i+1}^k-2w_{i}^k+w_{i-1}^k}{\Delta x^2} \]
In Eq. (47) we do not specify the total number of grid points in \(t\) (as we did with the grid in \(x\)); \(t\) is an initial value variable since it is typically time, and is defined over the semiinfinite interval \(0 \leq t \leq \infty\ .\)
Note that Eq. (49) explicitly gives the solution at the advanced point in \(t\) (at \(k+1\)) and therefore it is an explicit finite difference approximation of Eq. (41). We can now consider using Eq. (49) to step forward from an initial condition (IC) required by Eq. (41). Here we take as the initial condition \[\tag{50} w(x,t=0)=Ae^{-(\mu/D)x}+B \]
where \(A, B, \mu\) are constants to be specified. The finite difference form of Eq. (50) is \[\tag{51} w_i^0=Ae^{-(\mu/D)x_i}+B \]
(note that \(k=0\) at \(t=0\)).
We must also specify two boundary conditions (BCs) for Eq. (41) (since it is second order in \(x\)). We will use the Dirichlet BC at \(x=0\) \[\tag{52} w(x=0,t)=Ae^{-(\mu/D)(-\mu t)}+B \]
for which the finite difference form is (note that \(i=1\) at \(x=0\)) \[\tag{53} w_1^k=Ae^{-(\mu/D)(-\mu t^k)}+B \]
We use the Neumann BC at \(x=1\)
\[\tag{54}
\dfrac{\partial w(x=1,t)}{\partial x} = A (-\mu/D) e^{-(\mu/D)(1-\mu t)}
\]
for which the finite difference form is (note that \(i=N\) at \(x=1\)) \[\tag{55} w_{N+1}^k=w_{N-1}^k+2\Delta xAe^{-(\mu/D)(-\mu t^k)}+B \]
where \(w_{N+1}\) is a fictitious value that is outside the interval \(0 \leq x \leq 1\ ;\) it can be used to eliminate \(w_{N}\) in Eq. (49) for \(i=N\ .\)
Equations (49), (51), (53) and (55) constitute the full system of equations for the calculation of the numerical solution to Eq. (41). Note that we have replaced the original PDE, Eq. (41), with a set of approximating algebraic equations (Eqs. (49), (51), (53) and (55)) which can easily be programmed for a computer.
Also, an analytical solution to Eq. (41) (see particular solutions to the heat equation) can be used to evaluate the numerical solution \[\tag{56} w(x,t)=Ae^{(1/D)(\mu^2 t \pm \mu x)}+B \]
Equation (56) can be stated in the alternative form \[\tag{57} w(x,t)=Ae^{-(\mu/D)(x-\mu t)}+B \]
which corresponds to a traveling wave solution since \(x\) and \(t\) appear in the combination \(x-\mu t\ .\)
A short MATLAB program is listed in Appendix 1 based on Eqs. (49), (51), (53) and (55). Representative output from this program that compares the numerical solution from Eqs. (49), (51), (53) and (55) with the analytical solution, Eq. (57), indicates that the two solutions are in agreement to five figures, as reflected in Table 1.
| t = 0.05 w(x=xl/2,t) = 1.6376 wa(x=xl/2,t) = 1.6376 |
| t = 0.10 w(x=xl/2,t) = 1.6703 wa(x=xl/2,t) = 1.6703 |
| t = 0.15 w(x=xl/2,t) = 1.7047 wa(x=xl/2,t) = 1.7047 |
| t = 0.20 w(x=xl/2,t) = 1.7408 wa(x=xl/2,t) = 1.7408 |
| t = 0.25 w(x=xl/2,t) = 1.7788 wa(x=xl/2,t) = 1.7788 |
| t = 0.30 w(x=xl/2,t) = 1.8187 wa(x=xl/2,t) = 1.8187 |
| t = 0.35 w(x=xl/2,t) = 1.8607 wa(x=xl/2,t) = 1.8607 |
| t = 0.40 w(x=xl/2,t) = 1.9048 wa(x=xl/2,t) = 1.9048 |
| t = 0.45 w(x=xl/2,t) = 1.9512 wa(x=xl/2,t) = 1.9512 |
| t = 0.50 w(x=xl/2,t) = 2.0000 wa(x=xl/2,t) = 2.0000 |
The parameters that produced the numerical output in Table 1 are listed in Table 2.
| Parameter | Description | Value |
| \(D\) | diffusivity in Eq. (41) | 1 |
| \(\mu\), \(A\), \(B\) | constants in Eqs. (50)-(57) | 1, 1, 1 |
| \(N\), \(x_l\) | number of grid points and length in \(x\) | 21, 1 |
| \(h\), \(t_f\) | grid spacing in \(t\), final value of \(t\) | 0.001, 0.5 |
Additional parameters follow from the values in Table 2. Thus, the grid spacing in \(x\) is \(\Delta x=1/(21-1)=0.05\ .\) The number of steps in \(t\) taken along the solution is \(0.5/0.001=500\ .\)
Some of the parameters, particularly \(N\) and \(h\), were determined by trial and error to achieve a numerical solution of acceptable accuracy (e.g., five significant figures). We can note two additional points about these values:
- Acceptable values of \(N\) and \(h\) could be determined by observing the errors in the numerical solution through comparison with the exact solution (Eq. (57)) as illustrated in Table 1. However, in most PDE applications, an analytical solution is not available for assessing the accuracy of the numerical solution, and in fact, the motivation for using a numerical method is generally to produce a solution when an analytical solution is not available. In this case (no analytical solution), a useful procedure for estimating the numerical accuracy is to compute solutions for two different values of \(N\) and compare the numerical values. If the two solutions do not agree to an acceptable level, a third solution is computed with a still larger \(N\) (smaller grid spacing in \(x\)) and again the solutions are compared. Eventually, if spatial convergence is achieved, successive solutions will agree to the required accuracy. In this case, the accuracy of the numerical solution is inferred, but the exact error is not computed (or even known) when an analytical solution is not available. The same reasoning can be applied for temporal convergence with respect to \(t\), i.e., \(h\) is reduced until the successive solutions agree to a specified level.
- The preceding discussion was directed to achieving acceptable accuracy. Additionally, the values of \(N\) and \(h\) were selected to achieve a stable solution. The criterion for stability in the case of Eq. (49) (for parabolic PDE (41)) is \(\alpha = D\Delta t/\Delta x^2<1/2\ .\) For the solution of Table 1, \(\alpha=1\times0.001/0.05^2=0.4<0.5\ .\) If the critical value of the dimensionless parameter \(\alpha=0.5\) had been exceeded, the numerical solution would have become unstable (as manifest in numerical values of ever increasing magnitude). This conclusion can easily be confirmed by running the program of Appendix 1 for a choice of \(N\) and \(h\) for which \(\alpha>0.5\ .\) The stability constraint \(\alpha<0.5\) is a distinctive feature of the explicit finite difference approximation of Eq. (49). Thus, as \(\Delta x\) is reduced (\(N\) increased) to achieve better accuracy in the numerical solution, \(h\) must also be reduced to maintain stability (\(\alpha=Dh/\Delta x^2<0.5\)).
Finally, some descriptive comments about the details of the program in Appendix 1 are given immediately after the program listing.
Hyperbolic PDE
We now consider a numerical solution to the one-dimensional hyperbolic wave equation \[\tag{58} \frac{\partial^2w}{\partial t^2}-\frac{\partial^2w}{\partial x^2}=0 \]
We again have an analytical solution to evaluate the numerical solution. First, we include a velocity, \(c\), in the equation: \[\tag{59} \dfrac{\partial^2w}{\partial t^2} = c^2\dfrac{\partial^2w}{\partial x^2} \]
Note that \(c\) has the MKS units of m/s as expected and as inferred from Eq. (59) (note the units of the derivatives in \(x\) and \(t\)).
Since Eq. (59) is second order in \(x\) and \(t\), it requires two ICs and two BCs. We will take these as: \[\tag{60} w(x,t=0)=e^{-\lambda x^2} \]
\[\tag{61} \dfrac{\partial w(x,\,t=0)}{\partial t}=0 \]
\[\tag{62} w(x\to\infty,\,t)=0 \]
\[\tag{63} w(x\to-\infty,\,t)=0 \]
Equation (60) indicates the IC is a Gaussian pulse with the positive constant \(\lambda\) to be specified. Equation (61) indicates \(w(x,t)\) starts with zero "velocity". Equations (62)–(63) indicate that the solution \(w(x,t)\) does not depart from the initial value of zero specified by IC (60)–(61). In other words, \(\lambda\) is chosen large enough that the IC is effectively zero at \(x=\pm\infty\) and remains at this value for subsequent \(t\ .\)
An important difference between the parabolic problem of Eq. (41) and the hyperbolic problem of Eq. (59) is that the former is first order in \(t\) while the latter is second order in \(t\ .\) Therefore, in order to develop a numerical method for Eq. (58) (or Eq. (59)), we need an algorithm that can accommodate second derivatives in \(t\ .\) While such algorithms do exist, they generally are not required. Rather, we can express a PDE second order in \(t\) as two PDEs first order in \(t\ .\) For example, Eq. (59) can be written as \[\tag{64} \dfrac{\partial w}{\partial t} = u \]
\[\tag{65} \dfrac{\partial u}{\partial t} = c^2\dfrac{\partial^2w}{\partial x^2} \]
Equations (64)–(65) are first order in \(t\), and therefore an integration algorithm for first order equations, such as the Euler method of Eq. (48), can be used to move \(w(x,t)\) and \(u(x,t)\) forward in \(t\ .\) Thus, the fully discretized form of Eqs. (64)–(65) can be written as \[\tag{66} w_i^{k+1}=w_i^k+hu_{i}^k \]
\[\tag{67} u_{i}^{k+1}=u_{i}^k+hc^2\dfrac{w_{i+1}^k-2w_{i}^k+w_{i-1}^k}{\Delta x^2} \]
ICs (60)–(61) become (with \(k=0\) corresponding to \(t=0\))
\[\tag{68}
w_i^0=e^{-\lambda x_i^2}
\]
\[\tag{69} u_i^0=0 \]
For BCs (62)–(63), the infinite interval \(-\infty\leq x\leq\infty\) must be replaced by a finite
one \(-x_l\leq x\leq x_l\) (since computers can accommodate only finite numbers) where \(x_l\)
is selected so that it is effectively infinite; that is, the solution \(w(x,t)\) does not depart from
IC (60) at \(x=-x_l,x_l\) for \(t>0\ .\) The value of \(x_l\) and the corresponding number of grid points
in \(x\), \(N\), are specified subsequently.
The finite difference approximations of BCs (62)–(63) are \[\tag{70} w_1^k=0 \]
\[\tag{71} w_N^k=0 \]
Equations (66), (67), (68), (69), (70), (71) constitute the complete finite difference approximation of Eqs. (64), (65), (60), (61), (62), (63). A small MATLAB program for this approximation is given in Appendix 2, including some descriptive comments immediately after the program.
The general analytical solution to Eq. (59) can be written as (see also here) \[\tag{72} w(x,t)={\textstyle\frac12}[\varphi(x+ct)+\psi(x-ct)] \]
Let us take \(\varphi = \psi\) in the form of the Gaussian pulse of Eq. (60), i.e., \[\tag{73} w(x,t)={\textstyle\frac12}\bigl[e^{-\lambda(x+ct)^2}+e^{-\lambda(x-ct)^2}\bigr] \]
The plotted output from the program is given in Figure 1 and includes both the numerical solution of Eqs. (66), (67), (68), (69), (70), (71) and the analytical solution of Eq. (73).
We can note the following points about Figure 1:
- The initial Gaussian pulse at \(t=0\) (centered at \(x=0\) with unit maximum value) splits into two pulses traveling left and right with velocity \(c=1\) and maximum value of 0.5 according to Eq. (73).
- The pulses traveling left are centered at \(x=-10,-20,-30\) corresponding to \(t=10,20,30\) since \(c=1\ .\) The pulses traveling right are centered at \(x=10,20,30\) corresponding to \(t=10,20,30\ .\) These properties are characteristic of the traveling wave functions of Eq. (73) with arguments \(x+ct\) and \(x-ct\), respectively. In fact, the use of the word characteristic is particularly appropriate since the relations \(x+ct=\kappa\) and \(x-ct=\kappa\), where \(\kappa\) is a constant, are termed the characteristics of Eq. (59). Note that at the points along the solution given by these characteristics, the solution has a constant value. For example, for the peak values in Figure 1, \(\kappa=0\) and the solution is constant at the peak value 0.5 for \(x+ct=x-ct=0\ .\)
- The solution remains at zero for \(x=x_l=50\) so that the interval \(-50\leq x\leq50\) is equivalent to the infinite interval \(-\infty\leq x\leq\infty\ .\)
The parameters that produced the numerical output in Figure 1 are listed in Table 3.
| Parameter | Description | Value |
| \(c\) | velocity in Eq. (59) | 1 |
| \(\lambda\) | constant in Eqs. (60), (68), (73) | 0.05 |
| \(N\), \(x_l\) | number of grid points and half length in \(x\) | 201, 50 |
| \(h\), \(t_f\) | grid spacing in \(t\), final value of \(t\) | 0.0025, 30 |
Additional parameters follow from the values in Table 3. Thus, the grid spacing in \(x\) is \(\Delta x=2\times50/(201-1)=0.5\ .\) The number of steps in \(t\) taken along the solution is \(30/0.0025=12000\), which is large, but was selected to achieve good accuracy in the numerical solution.
To explore the accuracy of the numerical solution, we can consider the peak values of \(w(x,t)\) in Figure 1. This is a stringent test of the numerical solution since the curvature of the solution is greatest at these peaks. The results are summarized in Table 4.
| \(t\) | \(-x,\ x\) | Peak values |
| 0 | 0 | 1 |
| 10 | -10, 10 | 0.5006, 0.5006 |
| 20 | -20, 20 | 0.5011, 0.5011 |
| 30 | -30, 30 | 0.5015, 0.5015 |
Of course, the peak analytical values given by Eq. (73) are \(w(x=0,\,t=0)=1\ ,\) \(w(x=\pm t,\,t>0)=0.5\ .\) Table 4 indicates these peak values were attained within 0.5015 even for the largest value of \(t\ (=30)\) so that errors did not accumulate excessively as the solution progressed through \(t\ .\) These errors could be reduced by increasing the number of grid points in \(x\) above \(N=201\ .\) These errors could also presumably be reduced by using a more accurate (higher order) finite difference equation than the second order approximation of Eq. (45); to this end, MATLAB routines for second, fourth, sixth, eighth and tenth finite difference approximations are given in Hamdi, et al.
This explicit finite difference numerical solution also has a stability limit (like the preceding parabolic problem). In this case \[\tag{74} c\Delta t/\!\Delta x<1 \]
Stability constraint (74) is the Courant–Friedrichs–Lewy (or CFL) condition. For the present numerical solution, \(1\times0.0025/0.5=0.005\) so the CFL condition is easily satisfied. In other words, the parameters of Table 3 were chosen primarily for accuracy and not stability.
Next, we consider an elliptic problem.
Elliptic PDE
The elliptic PDE (Laplace's equation) \[\tag{75} \dfrac{\partial^2 w}{\partial x^2}+\dfrac{\partial^2 w}{\partial y^2}=0 \]
has two boundary value independent variables, \(x\) and \(y\), and no initial value variable. Thus, since the preceding numerical methods required an integration with respect to an initial value variable (\(t\)), and was accomplished by Euler's method, Eq. (48), we cannot develop a numerical solution for Eq. (75) directly using these methods. Rather, we will convert Eq. (75) from an elliptic problem to an associated parabolic problem by adding a derivative in an initial value variable. Thus, the PDE we will now consider is \[\tag{76} \dfrac{\partial w}{\partial t}=\dfrac{\partial^2 w}{\partial x^2}+\dfrac{\partial^2 w} {\partial y^2} \]
The idea then is to integrate Eq. (76) forward in \(t\) until the numerical solution approaches the condition \(\frac{\partial w}{\partial t}\to 0\) so that Eq. (76) then reverts back to Eq. (75), i.e., the solution under this condition is for Eq. (75) as required.
Equation (76) is first order in \(t\), and second order in \(x\) and \(y\ .\) It therefore requires one IC and two BCs (for \(x\) and \(y\)). For the IC, since \(t\) has been added to the problem only to provide a solution to Eq. (75) through Eq. (76), the choice of an IC is completely arbitrary (it is not part of the original problem). For the present analysis, we will use \[\tag{77} w(x,\,y,\,t=0)=\kappa \]
where \(\kappa\) is a constant to be selected (logically, it should be in the neighborhood of the expected solution to Eq. (75), but it's precise value is not critical to the success of the numerical method).
For the two BCs in \(x\), we will use homogeneous (zero) Dirichlet BCs \[\tag{78} w(x=0,\,y,\,t)=0 \]
\[\tag{79} w(x=x_l,\,y,\,t)=0 \]
where \(x_l\) is the upper boundary value of \(x\ .\)
For the first BC in \(y\), we will use a homogeneous Neumann BC \[\tag{80} \dfrac{\partial w(x,\,y=0,\,t)}{\partial y}=0 \]
For the second BC in \(y\), we will use a nonhomogeneous Neumann BC
\[\tag{81}
\dfrac{\partial w(x,\,y=y_l,\,t)}{\partial y} = \sin (\pi x)\,\pi \sinh (\pi y_l)
\]
where \(y_l\) is the upper boundary value of \(y\ .\) Note that \(\dfrac{\partial w(x,y=y_l)} {\partial y}\) in BC (81) is a function of \(x\ .\)
The analytical solution to Eqs. (75) (note, not Eq. (76)), (78)–(81) is a special case of one of the solutions stated previously \[\tag{82} w(x,y) = \sin (\pi x)\cosh(\pi y) \]
The analytical solution of Eq. (82) will be used to evaluate the numerical solution.
To develop a numerical solution to Eq. (76), we first replace all of the derivatives with finite differences in analogy with the preceding numerical solutions. The positions in \(x\), \(y\) and \(t\) will be denoted with indices \(i\), \(j\) and \(k\), respectively. The corresponding increments are \(\Delta x\ ,\) \(\Delta y\) and \(h\ .\) Application of these ideas to Eq. (76) gives the finite difference approximation \[ w_{i,j}^{k+1}=w_{i,j}^{k}+h\left (\dfrac{w_{i+1,j}^k-2w_{i,j}^k+w_{i-1,j}^k}{\Delta x^2} +\dfrac{w_{i,j+1}^k-2w_{i,j}^k+w_{i,j-1}^k}{\Delta y^2} \right ) \] \[\tag{83} i=1,2,\dots,N_x;\quad j=1,2,\dots,N_y;\quad k=1,2,\dots \]
BCs (78)–(79) in the difference notation are
\[\tag{84}
w_{1,j}^k = 0
\]
\[\tag{85} w_{N_x,j}^k = 0 \]
BC (80) in difference notation is
\[\tag{86}
w_{i,0}^k = w_{i,2}^k
\]
where \( w_{i,0}^k\) is a fictitious value that can be used in Eq. (83) for \(j=1\ .\)
BC (81) in difference notation is \[\tag{87} w_{i,N_y+1}^k = w_{i,N_y-1}^k+2\,\Delta y\sin (\pi x_i)\,\pi \sinh (\pi y_{N_y}) \]
where \( w_{i,N_y+1}^k\) is a fictitious value that can be used in Eq. (83) for \(j=N_y\ .\)
A short MATLAB program for Eqs. (83), (84), (85), (86) and (87) is listed in Appendix 3, followed by some explanatory comments. The output from this program is listed in Table 5.
| \(t\) | \(w(x_l/2,\,y_l/2,\,t)\) |
| 0.10 | 1.6254 |
| 0.20 | 2.1785 |
| 0.30 | 2.3883 |
| 0.40 | 2.4664 |
| 0.50 | 2.4956 |
| 0.60 | 2.5064 |
| 0.70 | 2.5105 |
| 0.80 | 2.5120 |
| 0.90 | 2.5125 |
| 1.00 | 2.5127 |
The midpoint values \(w(x=x_l/2,\,y=y_l/2,\,t)\) are listed in Table 5. The convergence of the solution of Eq. (76) to that of Eq. (75) is apparent, e.g., at \(t=1\ ,\) \(w(x=x_l/2,\,y=y_l/2,\,t)=2.5127\) while the value from the analytical solution, Eq. (82), is \(w_a(x=x_l/2,\,y=y_l/2)=2.5092\ ;\) the difference is \((2.5127-2.5092)/2.5092 \times 100 =0.14\%\ .\) Of course, the agreement will not be perfect because of the approximations used in Eqs. (83), (84), (85), (86) and (87). The parameters that produced the numerical output in Table 5 are listed in Table 6.
| Parameter | Description | Value |
| \(N_x\), \(x_l\) | number of grid points and half length in \(x\) | 21, 1 |
| \(N_y\), \(y_l\) | number of grid points and half length in \(y\) | 21, 1 |
| \(h\), \(t_f\) | grid spacing in \(t\), final value of \(t\) | 0.0005, 1 |
Additional parameters follow from the values in Table 6. Thus, the grid spacing in \(x\) is \(\Delta x = 1/(21-1) = 0.05\ .\) Similarly the spacing in \(y\) is \(\Delta y = 0.05\ .\) The spacing in \(x\) and \(y\) could be different if the variation of the solution \(u(x,y,t)\) in one direction is substantially greater than in the other; in other words, some tuning of the 2D grid in \(x\) and \(y\) could be required. The number of steps in \(t\) taken along the solution is \(1/0.0005 = 2000\), which is large, but was selected to achieve good accuracy in the numerical solution.
The stability constraint for the 2D problem of Eq. (76) (with a constant \(D\) multiplying the derivatives in \(x\) and \(y\)) is \[ D\left (\dfrac{1}{\Delta x^2}+\dfrac{1}{\Delta y^2}\right )h < \dfrac12 \] For the preceding solution, \(1\times(1/0.05^2 + 1/0.05^2)\times0.0005=0.4\) which meets the stability constraint.
The procedure of adding a derivative in \(t\) to elliptic Eq. (75), then solving the resulting parabolic problem, Eq. (76), is termed the method of false transients or the method of pseudo transients to suggest that the parabolic problem appears to have a transient phase that is not part of the original elliptic problem. The actual path that the parabolic problem takes to the solution of the elliptic problem is not relevant so long as the parabolic solution converges to the elliptic solution. In the present case, \(\kappa\) in IC (77) was taken as 1. Some other trial values indicated that this initial value is not critical (but it should be as close to the final value, for example 2.5092, as knowledge about the final value can provide).
Also, \(t\) can be considered a parameter in the solution of Eq. (76). This parametrization is an example of continuation in which the solution is continued from the given (assumed) starting value of Eq. (77) to the final value of 2.5127. The concept of continuation has been applied in many forms (and not just through the addition of a derivative as in Eq. (76)).
In general, the errors in the numerical solution of PDEs can result from the limited accuracy of all of the approximations used in the calculation. For example, the 0.14% error in the preceding solution can result from the discretization errors in \(x\), \(y\) and \(t\) in Eqs. (83), (84), (85), (86) and (87). In formulating a numerical method or algorithm for the solution of a PDE problem, it is necessary to balance the discretization errors so that one source of error does not dominate, and generally degrade, the numerical solution. This might, for example, require the choice of balanced values for \(N_x\), \(N_y\) and \(h\ .\) The effect of the numbers of discretization points generally can be inferred by varying these numbers and observing the effect on the numerical solution.
Remarks
Thus, control of approximation errors is central to the calculation of a numerical solution of acceptable accuracy. In the preceding examples, this control of errors can be accomplished in three ways:
- The number of grid points can be varied. Since in the numerical analysis literature, the grid spacing is often given the symbol \(h\) (as in Eq. (48)), systematic variation of the grid spacing to investigate accuracy is usually termed h refinement.
- The order of the approximation can be varied. For example, Euler's method of Eq. (48) is \(O(h^2)\) (second order correct) for one step or \(O(h)\) (first order correct) for a series of steps over a complete solution. In general, if the approximation is of order \(O(h^p)\), it is termed \(p\)th order. The symbol \(p\) is frequently used in the numerical analysis literature, and varying the order \(p\) to investigate accuracy is referred to as p refinement.
- The discretization errors can possibly be estimated and where the error is considered too large, the numerical algorithm can automatically insert grid points. Beyond that, the grid spacing does not have to be uniform, e.g., \(h\) in Eq. (48) does not have to remain constant throughout the numerical solution, and the numerical algorithm can automatically vary the spacing to concentrate the grid points where they are needed to achieve acceptable accuracy. This form of refinement is termed usually r refinement; we mention r refinement only to initiate possible interest in more advanced numerical methods.
The three preceding numerical solutions were developed using basic finite differences such as in Eqs. (45) and (48). However, many approaches to approximating derivatives in PDEs have been developed and used. Among these are finite elements, finite volumes, weighted residuals, e.g., collocation, Galerkin and spectral methods. Each of these methods has advantages and disadvantages, often according to the characteristics of the problem of interest (starting with the parabolic, hyperbolic and elliptic geometric classifications). Thus, an extensive literature for the numerical solution of PDEs is available, and we have only presented here a few basic concepts and examples.
The principal advantage of numerical methods applied to PDEs is that, in principle, PDEs of any number and complexity can be solved which is particularly useful when analytical solutions are not available. To illustrate how this might be done, the heat conduction equation with a nonlinear source term \(f(w)\) can easily be programmed by extending Eq. (49) to \[\tag{88} w_i^{k+1}=w_i^k+h\left [D\dfrac{w_{i+1}^k-2w_{i}^k+w_{i-1}^k}{\Delta x^2}+f(w_i^k)\right ] \]
where \(f(w)\) could be a MATLAB function to compute the nonlinearity in \(w(x,t)\ .\) Alternatively, the nonlinearity could be programmed directly in Eq. (88). For example, for a second order reaction, \[\tag{89} w_i^{k+1}=w_i^k+h\left[ D\dfrac{w_{i+1}^k-2w_{i}^k+w_{i-1}^k}{\Delta x^2}+(w_i^k)^2\right] \]
Note in particular how easily the nonlinearity is programmed.
As another example, a solution to the Burgers equation could be computed by extending Eq. (49) (with \(D = 1\)) to \[\tag{90} w_i^{k+1}=w_i^k+h\left [\dfrac{w_{i+1}^k-2w_{i}^k+w_{i-1}^k}{\Delta x^2}- w_i^k\dfrac{w_{i+1}^k-w_{i-1}^k}{2\,\Delta x}\right ] \]
where we have used the second order central finite difference approximation \[ \dfrac{\partial w}{\partial x}\approx\dfrac{w_{i+1}^k-w_{i-1}^k}{2\,\Delta x} \] Again, the nonlinear term in the Burgers equation is easily programmed as \[ w\dfrac{\partial w}{\partial x}\approx w_i^k\dfrac{w_{i+1}^k-w_{i-1}^k}{2\,\Delta x} \]
Finally, to conclude this introductory discussion of the numerical solution of PDEs, we used in the various examples integration with respect to an initial value variable, \(t\), by using the Euler method of Eq. (48). While Euler's method is general with respect to the form of the initial value integration, it does have two important limitations:
- The accuracy is first order (\(O(h)\)) and is therefore exact only for functions that vary linearly in \(t\ ;\) for higher order variations in \(t\), it is only approximate and generally requires a small \(h\) to achieve acceptable accuracy. One way around this accuracy limitation is to use a higher order integration algorithm in \(t\) such as the Runge-Kutta method (e.g., see the article by Shampine et al).
- Stability limits the step \(h\), such as in the stability constraint for Eq. (49), \(\alpha = D\Delta t/\Delta x^2 < 1/2\ .\) Such stability constraints for explicit integration methods such as Eq. (48) can be circumvented by using an implicit integration method, often termed a stiff method. We will not pursue various types of initial value integration methods here, but additional information is available (Shampine et al).
Thus, the Euler method is limited by both accuracy and stability. To circumvent these constraints, initial value integrators are used which are higher order (for improved accuracy) and stable (to provide a larger stable step \(h\) that may actually be unlimited, depending on the integrator and the application). As might be expected, such higher order methods are more complicated than the Euler method, but fortunately, they have been programmed in library routines that can easily be called and used. As an example, the MATLAB integrator ode15s has the following features:
- Good stability (it is a stiff integrator). The enhanced stability is achieved through implicit integration that requires the solution of algebraic equations at each step along the solution (each step of length \(h\)). An option for ode15s in the solution of the algebraic equations is the use of sparse matrix methods that are very efficient for large algebraic systems.
- Variable step and order, so that it automatically performs \(h\) and \(p\) refinement as the solution evolves.
- Library routines such as ode15s are written by experts who have included features that make them robust and reliable. Further, these routines have been used extensively and are therefore thoroughly tested.
The use of library routines for initial value integration is the basis for much of the work in the numerical method of lines solution of PDEs. This approach has been applied to a broad spectrum of PDE problems in 1D, 2D and 3D, including all of the major classes of PDEs (e.g., parabolic, hyperbolic and elliptic) in various orthogonal coordinate systems (e.g., Cartesian, cylindrical and spherical). The use of quality library routines provides an important step in the timely development of computer codes for new applications of PDEs whereby the analyst can take advantage of the work of experts, which is generally much more efficient and reliable than developing codes starting with just a general programming language.
Appendix 1. MATLAB program for a parabolic PDE
The MATLAB program that produced the results described in the Parabolic PDE for equation (41) follows.
clear all; clc
D=1; mu=1; A=1; B=1;
nx=21; xl=1; dx=0.05; x=[0:0.05:xl]; nx2=11;
nt=50; nout=10; h=0.001; t=0;
for i=1:nx; w(i)=A*exp(1/D*(mu^2*t-mu*x(i)))+B; end
for i1=1:nout
for i2=1:nt
for i=1:nx
if(i==1)w(1)=A*exp(1/D*(mu^2*t-mu*x(1)))+B; wt(1)=0;
elseif(i==nx)wf=w(nx-1)+(2*dx)*A*(-mu/D)*exp(1/D*(mu^2*t-mu*x(nx)));
wt(nx)=D*(wf-2*w(nx)+w(nx-1))/dx^2;
else wt(i)=D*(w(i+1)-2*w(i)+w(i-1))/dx^2; end
end
w=w+wt*h; t=t+h;
end
fprintf('t = %4.2f w(x=xl/2,t) = %6.4f wa(x=xl/2,t) = %6.4f\n', ...
t, w(nx2), A*exp(1/D*(mu^2*t-mu*x(nx2)))+B);
end
Listing 1: MATLAB program for the numerical solution of Eqs. (41), (50), (52), and (54)
We can note the following points about this program:
- Any previous files are first cleared. Then the parameters for the numerical solution of Eq. (41) are defined numerically as in Table 2. Note that a grid in \(x\) is set up as \(x=0,\,0.05,\,0.10,\,\dots,\,1\ .\)
clear all; clc D=1; mu=1; A=1; B=1; nx=21; xl=1; dx=0.05; x=[0:0.05:xl]; nx2=11; nt=50; nout=10; h=0.001; t=0;
- IC (51) is then implemented in a for loop (note that t=0).
for i=1:nx; w(i)=A*exp(1/D*(mu^2*t-mu*x(i)))+B; end
- Three nested for loops step Eq. (49) through \(x\) and \(t\)
for i1=1:nout
for i2=1:nt
for i=1:nx
Specifically, the outer loop (in it) gives the nout = 10 outputs in Table 1. Within each output interval (0.05), 50 steps in \(t\) are taken by the loop in i2 (in this way, the 500 steps in \(t\) produce outputs at only 10 points as displayed in Table 1). Within the loop in i, \(x\) spans the interval \(0\leq x\leq 1\ .\)
- BC (53) is programmed (at i=1) as
if(i==1)w(1)=A*exp(1/D*(mu^2*t-mu*x(1)))+B; wt(1)=0;
Also, since \(w(x=0,t)\) is specified by BC (52), the \(t\) derivative of this boundary value is set to zero so that the integration in \(t\) at \(x=0\) does not move \(w(x=0,t)\) from its prescribed value, i.e., wt(1)=0;).
- BC (55) is used to calculate the fictitious point wf which is then used in the finite difference calculation of Eq. (46).
elseif(i==nx)wf=w(nx-1)+(2*dx)*A*(-mu/D)*exp(1/D*(mu^2*t-mu*x(nx)));
wt(nx)=D*(wf-2*w(nx)+w(nx-1))/dx^2;
- For the interior points in \(x\) (\(x \neq 0, 1\)), Eq. (46) is programmed.
else wt(i)=D*(w(i+1)-2*w(i)+w(i-1))/dx^2; end
end
The second end terminates the inner loop in i, so that all of the values of \(x\) are handled.
- Euler's method, Eq. (48), is then used to advance the solution through nt=50 steps of length \(h=0.001\) through the intermediate for loop in i2.
w=w+wt*h; t=t+h; end
Note the use of the matrix facility of MATLAB (w and wt are vectors of length nx = 21). The end completes the intermediate loop in i2.
- The numerical solution is displayed to produce the output in Table 1, including the analytical solution of Eq. (56)
fprintf('t = %4.2f w(x=xl/2,t) = %6.4f wa(x=xl/2,t) = %6.4f\n', ...
t, w(nx2), A*exp(1/D*(mu^2*t-mu*x(nx2)))+B);
end
The end terminates the outer loop in i1.
Appendix 2. MATLAB program for a hyperbolic PDE
The MATLAB program that produced the results described in the Hyperbolic PDE section for equation (58) (and Eq. (59)) follows.
clear all; clc
c=1;
nx=201; xl=50; dx=0.5; x=[-xl:dx:xl];
nt=4000; nout=3; h=0.0025; t=0;
for i=1:nx; w(i)=exp(-0.05*x(i)^2); wt(i)=0; ...
wplot(1,i)=w(i); waplot(1,i)=exp(-0.05*x(i)^2); end
for i1=1:nout
for i2=1:nt
for i=1:nx
if(i==1) w( 1)=0; wt( 1)=0; wtt( 1)=0;
elseif(i==nx)w(nx)=0; wt(nx)=0; wtt(nx)=0;
else wtt(i)=c^2*(w(i+1)-2*w(i)+w(i-1))/dx^2; end
end
w=w+wt*h; wt=wt+wtt*h; t=t+h;
end
for i=1:nx
wplot (i1+1,i)=w(i);
waplot(i1+1,i)=0.5*(exp(-0.05*(x(i)-t)^2)+exp(-0.05*(x(i)+t)^2));
fprintf('\n %6.2f %6.2f,%8.4f %8.4f',t, x(i),wplot(i1+1,i),waplot(i1+1,i));
end
end
plot(x,wplot,'-',x,waplot,'o'); xlabel('x'); ylabel('w(x,t)');
title('w(x,t), t=0,10,20,30, solid - num, o - anal')
Listing 2: MATLAB program for the numerical solution of Eqs. (59), (60)–(63)
We can note the following points about this program:
- Any previous files are first cleared. Then the parameters for the numerical solution of Eq. (58) (or Eq. (59)) are defined numerically as in Table 3. Note that a grid in \(x\) is set up as \(x=-50,-49.5,\dots,50\), a total of 201 points. The choice of the upper and lower limits for \(x\), was made so that domain in \(x\) is essentially infinite (it accommodates BCs (62)–(63), (70)–(71) as discussed subsequently).
clear all; clc c=1; nx=201; xl=50; dx=0.5; x=[-xl:dx:xl]; nt=4000; nout=3; h=0.0025; t=0;
for i=1:nx; w(i)=exp(-0.05*x(i)^2); wt(i)=0; ...
wplot(1,i)=w(i); waplot(1,i)=exp(-0.05*x(i)^2); end
Both \(w(x,t)\) and \(\dfrac{\partial w(x,t)}{\partial t}\) are programmed (as w and wt) which is required since Eq. (59) is second order in \(t\) (as explained subsequently). Also, the initial value of \(w(x,t=0)\), both numerical and analytical, is stored for subsequent plotting (so that the plot will include IC (60)).
for i1=1:nout
for i2=1:nt
for i=1:nx
Specifically, the outer loop (in it) gives the nout = 3 outputs in Table 3. Within each output interval (1), 4000 steps in \(t\) are taken by the loop in i2 (in this way, the 4000 steps in \(t\) produce outputs at only 3 points as displayed in Figure 1). Within the loop in i, \(x\) spans the interval \(-50\leq x\leq 50\ .\)
if(i==1) w( 1)=0; wt( 1)=0; wtt( 1)=0;
elseif(i==nx)w(nx)=0; wt(nx)=0; wtt(nx)=0;
Also, since \(w\) is specified by BCs (62)–(63), the first and second derivatives in \(t\) of these boundary values are set to zero so that the integration in \(t\) does not move \(w\) from its zero values.
- For the interior points in \(x\) (\(x \neq -50, 50\)), Eq. (65) (discretized in \(x\)) is programmed.
else wtt(i)=c^2*(w(i+1)-2*w(i)+w(i-1))/dx^2; end
end
The second end terminates the inner loop in i, so that all of the values of \(x\) are handled.
- Euler's method (Eq. (48)) is then used to advance the solution through nt=4000 steps of length \(h=0.0025\) through the intermediate for loop in i2 according to Eqs. (66) and (67).
w=w+wt*h; wt=wt+wtt*h; t=t+h; end
Note the use of the matrix facility of MATLAB (w, wt and wtt are vectors of length nx = 201). The end completes the intermediate loop in i2.
- The numerical and analytical solutions are then stored at \(t=10,20,30\) (three passes through the outer loop in i1) followed by plotting (of Figure 1).
for i=1:nx
wplot (i1+1,i)=w(i);
waplot(i1+1,i)=0.5*(exp(-0.05*(x(i)-t)^2)+exp(-0.05*(x(i)+t)^2));
fprintf('\n %6.2f %6.2f,%8.4f %8.4f',t, x(i),wplot(i1+1,i),waplot(i1+1,i));
end
end
plot(x,wplot,'-',x,waplot,'o'); xlabel('x'); ylabel('w(x,t)');
title('w(x,t), t=0,10,20,30, solid - num, o - anal')
The final end terminates the outer loop in i1. Also, as a word of caution, the fprintf statement displays the complete numerical solution at \(t=10,20,30\) so that the numerical solution can be examined in detail (the peak values reported in Table 4 were taken from this output); if the program in Listing 2 is executed, the fprintf statement could be converted to a comment to reduce the numerical output.
Appendix 3. MATLAB program for an elliptic PDE
The MATLAB program that produced the results described in the Elliptic PDE section for equation (75) follows.
clear all; clc
nx=21; xl=1; dx=0.05; x=[0:0.05:xl]; nx2=11; mu=pi/xl;
ny=21; yl=1; dy=0.05; y=[0:0.05:yl]; ny2=11;
nt=200; nout=10; h=0.0005; t=0;
w=ones(nx,ny);
for i1=1:nout
for i2=1:nt
for i=1:nx
for j=1:ny
if(i== 1) w( 1,j)=0; wt( 1,j)=0;
elseif(i==nx)w(nx,j)=0; wt(nx,j)=0;
elseif(j== 1)wt( i,1)=(w(i+1,j)-2*w(i,j)+w(i-1,j))/dx^2 ...
+2*(w(i,2)-w(i,1))/dy^2;
elseif(j==ny)wf=w(i,ny-1)+2*dy*mu*sin(mu*x(i))*sinh(mu*y(ny));
wt(i,ny)=(w(i+1,ny)-2*w(i,ny)+w(i-1,ny))/dx^2 ...
+(wf-2*w(i,ny)+w(i,ny-1))/dy^2;
else wt(i,j)=(w(i+1,j)-2*w(i,j)+w(i-1,j))/dx^2 ...
+(w(i,j+1)-2*w(i,j)+w(i,j-1))/dy^2;
end
end
end
w=w+wt*h;
t=t+h;
end
fprintf('t = %5.2f w(xl/2,yl/2,t) = %7.4f \n',t,w(nx2,ny2))
end
fprintf('\n wa(xl/2,yl/2) = %7.4f\n',sin(mu*x(nx2))*cosh(mu*y(ny2)));
Listing 3: MATLAB program for the numerical solution of Eqs. (76)–(81)
We can note the following points about this program:
- Any previous files are first cleared. Then the parameters for the numerical solution of Eq. (76) are defined numerically as in Table 6. Note grids in \(x\) and \(y\) are set up as \(x=0,0.05,0.10,\dots,1\), \(y=0,0.05,0.10,\dots,1\ .\)
clear all; clc nx=21; xl=1; dx=0.05; x=[0:0.05:xl]; nx2=11; mu=pi/xl; ny=21; yl=1; dy=0.05; y=[0:0.05:yl]; ny2=11; nt=200; nout=10; h=0.0005; t=0;
- IC (77) is then implemented as
w=ones(nx,ny);
\(k=1\) in Eq. (77) is programmed with the MATLAB utility ones over the entire grid in \(x\) and \(y\), a total of nx\(\,\times\,\)ny points.
- Four nested for loops step Eq. (83) through \(x\), \(y\) and \(t\)
for i1=1:nout
for i2=1:nt
for i=1:nx
for j=1:ny
Specifically, the outer loop (in it) gives the nout = 10 outputs in Table 5. Within each output interval (0.0005), 200 steps in \(t\) are taken by the loop in i2 (in this way, the 2000 steps in \(t\) produce outputs at only 10 points as displayed in Table 5). Within the loop in i, \(x\) spans the interval \(0\leq x\leq 1\) and within the loop in j, \(y\) spans the interval \(0\leq y\leq 1\ .\)
- BC (78) is programmed (at i=1 for j=1 to j=ny) as
if(i== 1) w( 1,j)=0; wt( 1,j)=0;
Also, since \(w(x=0,y,t)\) is specified by BC (78), the \(t\) derivative of this boundary value is set to zero so that the integration in \(t\) at \(x=0\) does not move \(w(x=0,y,t)\) from its prescribed value, i.e., wt(1,j)=0;).
- BC (79) is programmed (at i=nx for j=1 to j=ny) as
elseif(i==nx)w(nx,j)=0; wt(nx,j)=0;
- Equation (83) (with subsequent application of Euler's method) is programmed using BC (80) (at i=1 to i=nx for j=1) as
elseif(j== 1)wt( i,1)=(w(i+1,j)-2*w(i,j)+w(i-1,j))/dx^2 ...
+2*(w(i,2)-w(i,1))/dy^2;
Note that we have used a finite difference approximation of Eq. (80), w(i,0) = w(i,2), to eliminate the fictitious value w(i,0) in Eq. (83).
- Equation (83) (with subsequent application of Euler's method) is programmed using BC (81) (at i=1 to i=nx for j=ny) as
elseif(j==ny)wf=w(i,ny-1)+2*dy*mu*sin(mu*x(i))*sinh(mu*y(ny));
wt(i,ny)=(w(i+1,ny)-2*w(i,ny)+w(i-1,ny))/dx^2 ...
+(wf-2*w(i,ny)+w(i,ny-1))/dy^2;
Note that we have used a finite difference approximation of Eq. (81) to eliminate the fictitious value w(i,ny+1)=wf in Eq. (83).
- For the interior points in \(x\) (\(x \neq 0, 1\)) and \(y\) (\(x \neq 0, 1\)), Eq. (83) (with subsequent application of Euler's method) is programmed.
else wt(i,j)=(w(i+1,j)-2*w(i,j)+w(i-1,j))/dx^2 ...
+(w(i,j+1)-2*w(i,j)+w(i,j-1))/dy^2;
end
end
end
The second and third end statements terminate the inner loops in i and j, so that all of the values of \(x\) and \(y\) are handled.
- Euler's method (Eq. (48)) is then used to advance the solution through nt=200 steps of length \(h=0.0005\) through the intermediate for loop in i2.
w=w+wt*h;
t=t+h;
end
Note the use of the matrix facility of MATLAB (w and wt are arrays of dimensions nx = 21, ny = 21). The end completes the intermediate loop in i2.
- The numerical solution is displayed to produce the output in Table 5.
fprintf('t = %5.2f w(xl/2,yl/2,t) = %7.4f \n',t,w(nx2,ny2))
end
The end terminates the outer loop in i1.
- The analytical solution, Eq. (82), is then printed at the end (as reflected in Table 5).
fprintf('\n wa(xl/2,yl/2) = %7.4f\n',sin(mu*x(nx2))*cosh(mu*y(ny2)));
The programs in Listings 1, 2 and 3 have an essential feature that should be pointed out: only the operations of storing and retrieving numbers, arithmetic and looping are used (which are done naturally and very efficiently by computers). In other words, the integration of the three PDEs was not done by the usual analytical mathematics, but rather, was done essentially with arithmetic. This is the essence of the numerical method, that is, to replace mathematics of differentiation and integration with something much simpler, but performed naturally by computers with extraordinary speed and precision. This simplification was accomplished by replacing the partial derivatives in the PDEs with algebraic approximations, in this case, finite differences.
A central question then is how accurate the numerical solution is if approximations are used to compute it. We have observed the accuracy can be improved with more calculations (e.g., by reducing \(\Delta x\), \(\Delta y\) and \(h\)). In other words there is a tradeoff between accuracy and calculations, and we basically assume the computer can do enough calculations to achieve the required accuracy. Historically, the increasing speed of computers has produced solutions to PDE problems with acceptable accuracy and continually increasing complexity. Thus, for the future, we can expect to compute numerical solutions of PDEs with essentially unlimited complexity if sufficient computing power is available.
References
- R. Courant and D. Hilbert, Methods of Mathematical Physics. Volume 2. Partial Differential Equations, Wiley-VCH, 1989.
- L. C. Evans, Partial Differential Equations, American Mathematical Society, Providence, 1998.
- S. J. Farlow, Partial Differential Equations for Scientists and Engineers, Dover Publications Inc., 1993.
- F. John, Partial Differential Equations. Fourth Edition, Springer, 1991.
- J. Jost, Partial Differential Equations, Springer-Verlag, New York, 2002.
- I. G. Petrovskii, Partial Differential Equations, W. B. Saunders Co., Philadelphia, 1967.
- Y. Pinchover and J. Rubinstein, An Introduction to Partial Differential Equations, Cambridge University Press, Cambridge, 2005.
- A. D. Polyanin, Handbook of Linear Partial Differential Equations for Engineers and Scientists, Chapman & Hall/CRC Press, Boca Raton, 2002.
- A. D. Polyanin and V. F. Zaitsev, Handbook of Nonlinear Partial Differential Equations, Chapman & Hall/CRC Press, Boca Raton, 2004.
- A. D. Polyanin, V. F. Zaitsev, and A. Moussiaux, Handbook of First Order Partial Differential Equations, Taylor & Francis, London, 2002.
- D. L. Powers, Boundary Value Problems, Fifth Edition: and Partial Differential Equations, Elsevier Academic Press, 2005.
- W. E. Schiesser, Computational Mathematics in Engineering and Applied Science: ODEs, DAEs, and PDEs, CRC Press, Boca Raton, 1993.
- I. Stakgold, Boundary Value Problems of Mathematical Physics, Vols. I, II, SIAM, Philadelphia, 2000.
- A. N. Tikhonov and A. A. Samarskii, Equations of Mathematical Physics, Dover Publ., New York, 1990.
- D. Zwillinger, Handbook of Differential Equations (3rd edition), Academic Press, Boston, 1997.
External links
- Partial Differential Equations: Exact Solutions at EqWorld: The World of Mathematical Equations [1]
- Partial Differential Equations: Index of PDEs at EqWorld: The World of Mathematical Equations [2]
- Partial Differential Equations: Methods at EqWorld: The World of Mathematical Equations [3]
- Partial Differential Equation at Wolfram MathWorld by Eric Weisstein [4]
- Example problems with solutions at ExampleProblems.com [5]
- General reference for numerical methods at Scholarpedia [6]
- Introduction to numerical methods for partial differential equations at Scholarpedia [7]
