Calculus/Partial differential equations

First order

Any partial differential equation of the form

$h_{1}{\frac {\partial u}{\partial x_{1}}}+h_{2}{\frac {\partial u}{\partial x_{2}}}\cdots +h_{n}{\frac {\partial u}{\partial x_{n}}}=b$

where h₁, h₂ … h_n, and b are all functions of both u and Rⁿ can be reduced to a set of ordinary differential equations.

To see how to do this, we will first consider some simpler problems.

Special cases

We will start with the simple PDE

$u_{z}(x,y,z)=u(x,y,z)\quad (1)$

Because u is only differentiated with respect to z, for any fixed x and y we can treat this like the ODE, du/dz=u. The solution of that ODE is ce^z, where c is the value of u when z=0, for the fixed x and y

Therefore, the solution of the PDE is

$u(x,y,z)=u(x,y,0)e^{z}$

Instead of just having a constant of integration, we have an arbitrary function. This will be true for any PDE.

Notice the shape of the solution, an arbitrary function of points in the xy, plane, which is normal to the 'z' axis, and the solution of an ODE in the 'z' direction.

Now consider the slightly more complex PDE

$a_{x}u_{x}+a_{y}u_{y}+a_{z}u_{z}=h(u)\quad (2)$

where h can be any function, and each a is a real constant.

We recognize the left hand side as being a·∇, so this equation says that the differential of u in the a direction is h(u). Comparing this with the first equation suggests that the solution can be written as an arbitrary function on the plane normal to a combined with the solution of an ODE.

Remembering from Calculus/Vectors that any vector r can be split up into components parallel and perpendicular to a,

$\mathbf {r} =\mathbf {r} _{\perp }+\mathbf {r} _{\|}=\left(\mathbf {r} -{\frac {(\mathbf {r} \cdot \mathbf {a} )\mathbf {a} }{|\mathbf {a} |^{2}}}\right)+{\frac {(\mathbf {r} \cdot \mathbf {a} )\mathbf {a} }{|\mathbf {a} |^{2}}}$

we will use this to split the components of r in a way suggested by the analogy with (1).

Let's write

$\mathbf {r} =(x,y,z)=\mathbf {r} _{\perp }+s\mathbf {a} \quad s={\frac {\mathbf {r} \cdot \mathbf {a} }{\mathbf {a} \cdot \mathbf {a} }}$

and substitute this into (2), using the chain rule. Because we are only differentiating in the a direction, adding any function of the perpendicular vector to s will make no difference.

First we calculate grad s, for use in the chain rule,

$\nabla s={\frac {\mathbf {a} }{a^{2}}}$

On making the substitution into (2), we get,

$h(u)=\mathbf {a} \cdot \nabla s{\frac {d}{ds}}u(s)={\frac {\mathbf {a} \cdot \mathbf {a} }{\mathbf {a} \cdot \mathbf {a} }}{\frac {d}{ds}}u(s)={\frac {du}{ds}}$

which is an ordinary differential equation with the solution

$s=c(\mathbf {r} _{\perp })+\int ^{u}{\frac {dt}{h(t)}}$

The constant c can depend on the perpendicular components, but not upon the parallel coordinate. Replacing s with a monotonic scalar function of s multiplies the ODE by a function of s, which doesn't affect the solution.

Example:

$u(x,t)_{x}=u(x,t)_{t}$

For this equation, a is (1, -1), s=x-t, and the perpendicular vector is (x+t)(1, 1). The reduced ODE is du/ds=0 so the solution is

u=f(x+t)

To find f we need initial conditions on u. Are there any constraints on what initial conditions are suitable?

Consider, if we are given

• u(x,0), this is exactly f(x),
• u(3t,t), this is f(4t) and f(t) follows immediately
• u(t³+2t,t), this is f(t³+3t) and f(t) follows, on solving the cubic.
• u(-t,t), then this is f(0), so if the given function isn't constant we have a inconsistency, and if it is the solution isn't specified off the initial line.

Similarly, if we are given u on any curve which the lines x+t=c intersect only once, and to which they are not tangent, we can deduce f.

For any first order PDE with constant coefficients, the same will be true. We will have a set of lines, parallel to r=at, along which the solution is gained by integrating an ODE with initial conditions specified on some surface to which the lines aren't tangent.

If we look at how this works, we'll see we haven't actually used the constancy of a, so let's drop that assumption and look for a similar solution.

The important point was that the solution was of the form u=f(x(s),y(s)), where (x(s),y(s)) is the curve we integrated along -- a straight line in the previous case. We can add constant functions of integration to s without changing this form.

Consider a PDE,

$a(x,y)u_{x}+b(x,y)u_{y}=c(x,y,u)$

For the suggested solution, u=f(x(s),y(s)), the chain rule gives

${\frac {du}{ds}}={\frac {dx}{ds}}u_{x}+{\frac {dy}{ds}}u_{y}$

Comparing coefficients then gives

${\frac {dx}{ds}}=a(x,y)\quad {\frac {dy}{ds}}=b(x,y)\quad {\frac {du}{ds}}=c(x,y,u)$

so we've reduced our original PDE to a set of simultaneous ODE's. This procedure can be reversed.

The curves (x(s),y(s)) are called characteristics of the equation.

Example: Solve $yu_{x}=xu_{y}$ given u=f(x) for x≥0 The ODE's are

${\frac {dx}{ds}}=y\quad {\frac {dy}{ds}}=-x\quad {\frac {du}{ds}}=0$

subject to the initial conditions at s=0,

$x(0)=r\quad y(0)=0\quad u(0)=f(r)\quad r\geq 0$

This ODE is easily solved, giving

$x(s)=r\cos s\quad y(s)=r\sin s\quad u(s)=f(r)$

so the characteristics are concentric circles round the origin, and in polar coordinates u(r,θ)=f(r)

Considering the logic of this method, we see that the independence of a and b from u has not been used either, so that assumption too can be dropped, giving the general method for equations of this quasilinear form.

Quasilinear

Summarising the conclusions of the last section, to solve a PDE

$a_{1}(u,\mathbf {x} ){\frac {\partial u}{\partial x_{1}}}+a_{2}(u,\mathbf {x} ){\frac {\partial u}{\partial x_{2}}}\cdots +a_{n}(u,\mathbf {x} ){\frac {\partial u}{\partial x_{n}}}=b(u,\mathbf {x} )$

subject to the initial condition that on the surface, (x₁(r₁,…,r_n-1, …x_n(r₁,…,r_n-1), u=f(r₁,…,r_n-1) --this being an arbitrary paremetrisation of the initial surface--

• we transform the equation to the equivalent set of ODEs,

${\frac {dx_{1}}{ds}}=a_{1}\ \ldots \ {\frac {dx_{n}}{ds}}=a_{n}\quad {\frac {du}{ds}}=b$

subject to the initial conditions

$x_{i}(0)=f(r_{1},\ldots ,r_{n-1})\quad u=f(r_{1},r_{2},\ldots r_{n-1})$

• Solve the ODE's, giving x_i as a function of s and the r_i.
• Invert this to get s and the r_i as functions of the x_i.
• Substitute these inverse functions into the expression for u as a function of s and the r_i obtained in the second step.

Both the second and third steps may be troublesome.

The set of ODEs is generally non-linear and without analytical solution. It may even be easier to work with the PDE than with the ODEs.

In the third step, the r_i together with s form a coordinate system adapted for the PDE. We can only make the inversion at all if the Jacobian of the transformation to Cartesian coordinates is not zero,

${\begin{vmatrix}{\frac {\partial x_{1}}{\partial r_{1}}}&\cdots &{\frac {\partial x_{1}}{\partial r_{n-1}}}&a_{1}\\\vdots &\ddots &&\vdots \\{\frac {\partial x_{n}}{\partial r_{1}}}&\cdots &{\frac {\partial x_{n}}{\partial r_{n-1}}}&a_{n}\end{vmatrix}}\neq 0$

This is equivalent to saying that the vector (a₁, &hellip:, a_n) is never in the tangent plane to a surface of constant s.

If this condition is not false when s=0 it may become so as the equations are integrated. We will soon consider ways of dealing with the problems this can cause.

Even when it is technically possible to invert the algebraic equations it is obviously inconvenient to do so.

Example

To see how this works in practice, we will
a/ consider the PDE,

$uu_{x}+u_{y}+u_{t}=0$

with generic initial condition,

$u=f(x,y){\mbox{ on }}t=0$

Naming variables for future convenience, the corresponding ODE's are

${\frac {dx}{d\tau }}=u\quad {\frac {dy}{d\tau }}=1\quad {\frac {dt}{d\tau }}=1\quad {\frac {du}{d\tau }}=0\quad$

subject to the initial conditions at τ=0

$x=r\quad y=s\quad t=0\quad u=f(r,s)$

These ODE's are easily solved to give

$x=r+f(r,s)\tau \quad y=s+\tau \quad t=\tau \quad u=f(r,s)$

These are the parametric equations of a set of straight lines, the characteristics.

The determinant of the Jacobian of this coordinate transformation is

${\begin{vmatrix}1+\tau {\frac {\partial f}{\partial r}}&\tau {\frac {\partial f}{\partial s}}&f\\0&1&1\\0&0&1\end{vmatrix}}=1+\tau {\frac {\partial f}{\partial r}}$

This determinant is 1 when t=0, but if f_r is anywhere negative this determinant will eventually be zero, and this solution fails.

In this case, the failure is because the surface $sf_{r}=-1$ is an envelope of the characteristics.

For arbitrary f we can invert the transformation and obtain an implicit expression for u

$u=f(x-tu,y-x)$

If f is given this can be solved for u.

1/ $f(x,y)=ax$, The implicit solution is

$u=a(x-tu)\Rightarrow u={\frac {ax}{1+at}}$

This is a line in the u-x plane, rotating clockwise as t increases. If a is negative, this line eventually become vertical. If a is positive, this line tends towards u=0, and the solution is valid for all t.

For f(x,y)=x² The implicit solution is:

$u=(x-tu)^{2}\Rightarrow u={\frac {1+2tx-{\sqrt {1+4tx}}}{2t^{2}}}$

This solution clearly fails when $1+4tx<0$, which is just when $sf_{r}=-1$ . For any t > 0 this happens somewhere. As t increases this point of failure moves toward the origin.

Notice that the point where u=0 stays fixed. This is true for any solution of this equation, for all values of f

We will see later that we can find a solution after this time, if we consider discontinuous solutions. We can think of this as a shockwave.

3/ $f(x,y)=\sin(xy)$
The implicit solution is

$u(x,y,t)=\sin \left((x-tu)(y-x)\right)$

and we can not solve this explitely for u. The best we can manage is a numerical solution of this equation.

We can also consider the closely related PDE

$uu_{x}+u_{y}+u_{t}=y$

The corresponding ODE's are

${\frac {dx}{d\tau }}=u\quad {\frac {dy}{d\tau }}=1\quad {\frac {dz}{d\tau }}=1\quad {\frac {du}{d\tau }}=y\quad$

subject to the initial conditions at τ=0

$x=r\quad y=s\quad t=0\quad u=f(r,s)$

These ODE's are easily solved to give

$x=r+\tau f+{\frac {1}{2}}s\tau ^{2}+{\frac {1}{6}}\tau ^{3}\quad y=s+\tau \quad t=\tau \quad u=f+s\tau +{\frac {1}{2}}\tau ^{2}$

Writing f in terms of u, s, and τ, then substituting into the equation for x gives an implicit solution

$u(x,y,t)=f(x-ut+{\frac {1}{2}}yt^{2}-{\frac {1}{6}}t^{3},y-t)+yt-{\frac {1}{2}}t^{2}$

It is possible to solve this for u in some special cases, but in general we can only solve this equation numerically. However, we can learn much about the global properties of the solution from further analysis

Characteristic initial value problems

What if initial conditions are given on a characteristic, on an envelope of characteristics, on a surface with characteristic tangents at isolated points?

Discontinuous solutions

So far, we've only considered smooth solutions of the PDE, but this is too restrictive. We may encounter initial conditions which aren't smooth, e.g.

$u_{t}=cu_{x}\quad u(x,0)=\left\{{\begin{matrix}1,&x\geq 0\\0,&x<0\end{matrix}}\right.$

If we were to simply use the general solution of this equation for smooth initial conditions,

$u(x,t)=u(x+ct,0)\,$

we would get

$u(x,t)=\left\{{\begin{matrix}1,&x+ct\geq 0\\0,&x+ct<0\end{matrix}}\right.$

which appears to be a solution to the original equation. However, since the partial differentials are undefined on the characteristic x+ct=0, so it becomes unclear what it means to say that the equation is true at that point.

We need to investigate further, starting by considering the possible types of discontinuities.

If we look at the derivations above, we see we've never use any second or higher order derivatives so it doesn't matter if they aren't continuous, the results above will still apply.

The next simplest case is when the function is continuous, but the first derivative is not, e.g. |x|. We'll initially restrict ourselves to the two-dimensional case, u(x, t) for the generic equation.

$a(x,t)u_{x}+b(x,t)u_{t}=c(u,x,t)\quad (1)$

Typically, the discontinuity is not confined to a single point, but is shared by all points on some curve, (x₀(s), t₀(s) )

Then we have

${\begin{matrix}x>x_{0}&\lim _{x\rightarrow x_{0}}=u_{+}\\x<x_{0}&\lim _{x\rightarrow x_{0}}=u_{-}\end{matrix}}$

We can then compare u and its derivatives on both sides of this curve.

It will prove useful to name the jumps across the discontinuity. We say

$[u]=u_{+}-u_{-}\quad [u_{x}]=u_{x+}-u_{x-}\quad [u_{t}]=u_{t+}-u_{t-}$

Now, since the equation (1) is true on both sides of the discontinuity, we can see that both u₊ and u_-, being the limits of solutions, must themselves satisfy the equation. That is,

${\begin{matrix}a(x,t)u_{+x}+b(x,t)u_{+t}&=&c(u_{+},x,t)\\a(x,t)u_{-x}+b(x,t)u_{-t}&=&c(u_{-},x,t)\end{matrix}}{\mbox{ where }}{\begin{matrix}x=x_{0}(s)\\t=t_{0}(s)\end{matrix}}$

Subtracting then gives us an equation for the jumps in the differentials

$a(x,t)[u_{x}]+b(x,t)[u_{t}]=0\,$

We are considering the case where u itself is continuous so we know that [u]=0. Differentiating this with respect to s will give us a second equation in the differential jumps.

${\frac {dx_{0}}{ds}}[u_{x}]+{\frac {dt_{0}}{ds}}[u_{t}]=0$

The last two equations can only be both true if one is a multiple of the other, but multiplying s by a constant also multiplies the second equation by that same constant while leaving the curve of discontinuity unchanged, hence we can without loss of generality define s to be such that

${\frac {dx_{0}}{ds}}=a\quad {\frac {dt_{0}}{ds}}=b$

But these are the equations for a characteristic, i.e. discontinuities propagate along characteristics. We could use this property as an alternative definition of characteristics.

We can deal similarly with discontinuous functions by first writing the equation in conservation form, so called because conservation laws can always be written this way.

$(au)_{x}+(bu)_{t}=a_{x}u+b_{t}u+c\quad (1)$

Notice that the left hand side can be regarded as the divergence of (au, bu). Writing the equation this way allows us to use the theorems of vector calculus.

Consider a narrow strip with sides parallel to the discontinuity and width h

We can integrate both sides of (1) over R, giving

$\int _{R}(au)_{x}+(bu)_{t}\,dxdt=\int _{R}(a_{x}+b_{t})u+c\,dxdt$

Next we use Green's theorem to convert the left hand side into a line integral.

$\oint _{\partial R}audt-budx=\int _{R}(a_{x}+b_{t})u+c\,dxdt$

Now we let the width of the strip fall to zero. The right hand side also tends to zero but the left hand side reduces to the difference between two integrals along the part of the boundary of R parallel to the curve.

$\int au_{+}dt-bu_{+}dx-\int au_{-}dt-bu_{-}dx=0$

The integrals along the opposite sides of R have different signs because they are in opposite directions.

For the last equation to always be true, the integrand must always be zero, i.e.

$\left(a{\frac {dt_{0}}{ds}}-b{\frac {dx_{0}}{ds}}\right)[u]=0$

Since, by assumption [u] isn't zero, the other factor must be, which immediately implies the curve of discontinuity is a characteristic.

Once again, discontinuities propagate along characteristics.

Above, we only considered functions of two variables, but it is straightforward to extend this to functions of n variables.

The initial condition is given on an n-1 dimensional surface, which evolves along the characteristics. Typical discontinuities in the initial condition will lie on a n-2 dimensional surface embedded within the initial surface. This surface of discontinuity will propagate along the characteristics that pass through the initial discontinuity.

The jumps themselves obey ordinary differential equations, much as u itself does on a characteristic. In the two dimensional case, for u continuous but not smooth, a little algebra shows that

${\frac {d[u_{x}]}{ds}}=[u_{x}]\left({\frac {\partial c}{\partial u}}+a{\frac {b_{x}}{b}}-a_{x}\right)$

while u obeys the same equation as before,

${\frac {du}{ds}}=c$

We can integrate these equations to see how the discontinuity evolves as we move along the characteristic.

We may find that, for some future s, [u_x] passes through zero. At such points, the discontinuity has vanished, and we can treat the function as smooth at that characteristic from then on.

Conversely, we can expect that smooth functions may, under the right circumstances, become discontinuous.

To see how all this works in practice we'll consider the solutions of the equation

$u_{t}+uu_{x}=0\quad u(x,0)=f(x)$

for three different initial conditions.

The general solution, using the techniques outlined earlier, is

$u=f(x-tu)\,\!$

u is constant on the characteristics, which are straight lines with slope dependent on u.

First consider f such that

$f(x)=\left\{{\begin{matrix}1&x>a\\{\frac {x}{a}}&a\geq x>0\\0&x\leq 0\end{matrix}}\right.\quad a>0$

While u is continuous its derivative is discontinuous at x=0, where u=0, and at x=a, where u=1. The characteristics through these points divide the solution into three regions.

All the characteristics to the right of the characteristic through x=a, t=0 intersect the x-axis to the right of x=1, where u=1 so u is 1 on all those characteristics, i.e whenever x-t>a.

Similarly the characteristic through the origin is the line x=0, to the left of which u remains zero.

We could find the value of u at a point in between those two characteristics either by finding which intermediate characteristic it lies on and tracing it back to the initial line, or via the general solution.

Either way, we get

$f(x)=\left\{{\begin{matrix}1&x-t>a\\{\frac {x}{a+t}}&a+t\geq x>0\\0&x\leq 0\end{matrix}}\right.$

At larger t the solution u is more spread out than at t=0 but still the same shape.

We can also consider what happens when a tends to 0, so that u itself is discontinuous at x=0.

If we write the PDE in conservation form then use Green's theorem, as we did above for the linear case, we get

$[u]{\frac {dx_{0}}{ds}}={\frac {1}{2}}[u^{2}]{\frac {dt_{0}}{ds}}$

[u²] is the difference of two squares, so if we take s=t we get

${\frac {dx_{0}}{dt}}={\frac {1}{2}}(u_{-}+u_{+})$

In this case the discontinuity behaves as if the value of u on it were the average of the limiting values on either side.

However, there is a caveat.

Since the limiting value to the left is u_- the discontinuity must lie on that characteristic, and similarly for u₊; i.e the jump discontinuity must be on an intersection of characteristics, at a point where u would otherwise be multivalued.

For this PDE the characteristic can only intersect on the discontinuity if

$u_{-}>u_{+}\,$

If this is not true the discontinuity can not propagate. Something else must happen.

The limit a=0 is an example of a jump discontinuity for which this condition is false, so we can see what happens in such cases by studying it.

Taking the limit of the solution derived above gives

$f(x)=\left\{{\begin{matrix}1&x>t\\{\frac {x}{t}}&t\geq x>0\\0&x\leq 0\end{matrix}}\right.$

If we had taken the limit of any other sequence of initials conditions tending to the same limit we would have obtained a trivially equivalent result.

Looking at the characteristics of this solution, we see that at the jump discontinuity characteristics on which u takes every value between 0 and 1 all intersect.

At later times, there are two slope discontinuities, at x=0 and x=t, but no jump discontinuity.

This behaviour is typical in such cases. The jump discontinuity becomes a pair of slope discontinuities between which the solution takes all appropriate values.

Now, lets consider the same equation with the initial condition

$f(x)=\left\{{\begin{matrix}1&x\leq 0\\1-{\frac {x}{a}}&a\geq x>0\\0&x>a\end{matrix}}\right.\quad a>0$

This has slope discontinuities at x=0 and x=a, dividing the solution into three regions.

The boundaries between these regions are given by the characteristics through these initial points, namely the two lines

$x=t\quad x=a$

These characteristics intersect at t=a, so the nature of the solution must change then.

In between these two discontinuities, the characteristic through x=b at t=0 is clearly

$x=\left(1-{\frac {b}{a}}\right)t+b\quad 0\leq b\leq a$

All these characteristics intersect at the same point, (x,t)=(a,a).

We can use these characteristics, or the general solution, to write u for t<a

$u(x,t)=\left\{{\begin{matrix}1&x\leq t\\{\frac {a-x}{a-t}}&a\geq x>t\\0&x>a\end{matrix}}\right.\quad a>t\geq 0$

As t tends to a, this becomes a step function. Since u is greater to the left than the right of the discontinuity, it meets the condition for propagation deduced above, so for t>a u is a step function moving at the average speed of the two sides.

$u(x,t)=\left\{{\begin{matrix}1&x\leq {\frac {a+t}{2}}\\0&x>{\frac {a+t}{2}}\end{matrix}}\right.\quad t\geq a\geq 0$

This is the reverse of what we saw for the initial condition previously considered, two slope discontinuities merging into a step discontinuity rather than vice versa. Which actually happens depends entirely on the initial conditions. Indeed, examples could be given for which both processes happen.

In the two examples above, we started with a discontinuity and investigated how it evolved. It is also possible for solutions which are initially smooth to become discontinuous.

For example, we saw earlier for this particular PDE that the solution with the initial condition u=x² breaks down when 2xt+1=0. At these points the solution becomes discontinuous.

Typically, discontinuities in the solution of any partial differential equation, not merely ones of first order, arise when solutions break down in this way and propagate similarly, merging and splitting in the same fashion.

Fully non-linear PDEs

It is possible to extend the approach of the previous sections to reduce any equation of the form

$F(x_{1},x_{2},\ldots ,x_{n},u,u_{x_{1}},u_{x_{2}},\ldots ,u_{x_{n}})=0$

to a set of ODE's, for any function, F.

We will not prove this here, but the corresponding ODE's are

${\frac {dx_{i}}{d\tau }}={\frac {\partial F}{\partial u_{i}}}\quad {\frac {du_{i}}{d\tau }}=-\left({\frac {\partial F}{\partial x_{i}}}+u_{i}{\frac {\partial F}{\partial u}}\right)\quad {\frac {du}{d\tau }}=\sum _{i=1}^{n}u_{i}{\frac {\partial F}{\partial u_{i}}}$

If u is given on a surface parameterized by r₁…r_n then we have, as before, n initial conditions on the n, x_i

$\tau =0\quad x_{i}=f_{i}(r_{1},r_{2},\ldots ,r_{n-1})$

given by the parameterization and one initial condition on u itself,

$\tau =0\quad u=f(r_{1},r_{2},\ldots ,r_{n-1})$

but, because we have an extra n ODEs for the u_i's, we need an extra n initial conditions.

These are, n-1 consistency conditions,

$\tau =0\quad {\frac {\partial f}{\partial r_{i}}}=\sum _{j=1}^{n-1}u_{i}{\frac {\partial f_{i}}{\partial r_{j}}}$

which state that the u_i's are the partial derivatives of u on the initial surface, and one initial condition

$\tau =0\quad F(x_{1},x_{2},\ldots ,x_{n},u,u_{1},u_{2},\ldots ,u_{n})=0$

stating that the PDE itself holds on the initial surface.

These n initial conditions for the u_i will be a set of algebraic equations, which may have multiple solutions. Each solution will give a different solution of the PDE.

Example

Consider

$u_{t}=u_{x}^{2}+u_{y}^{2},\quad u(x,y,0)=x^{2}+y^{2}$

The initial conditions at τ=0 are

${\begin{matrix}x=r&y=s&t=0&u=r^{2}+s^{2}\\u_{x}=2r&u_{y}=2s&u_{t}=4(r^{2}+s^{2})&\end{matrix}}$

and the ODE's are

${\begin{matrix}{\frac {dx}{d\tau }}=-2u_{x}&{\frac {dy}{d\tau }}=-2u_{y}&{\frac {dt}{d\tau }}=1&{\frac {du}{d\tau }}=u_{t}-2(u_{x}^{2}+u_{y}^{2})\\{\frac {du_{x}}{d\tau }}=0&{\frac {du_{y}}{d\tau }}=0&{\frac {du_{t}}{d\tau }}=0&\end{matrix}}$

Note that the partial derivatives are constant on the characteristics. This always happen when the PDE contains only partial derivatives, simplifying the procedure.

These equations are readily solved to give

$x=r(1-4\tau )\quad y=s(1-4\tau )\quad t=\tau \quad u=(r^{2}+s^{2})(1-4\tau )$

On eliminating the parameters we get the solution,

$u={\frac {x^{2}+y^{2}}{1-4t}}$

which can easily be checked. abc

Second order

Suppose we are given a second order linear PDE to solve

$a(x,y)u_{xx}+b(x,y)u_{xy}+c(x,y)u_{yy}=d(x,y)u_{x}+e(x,y)u_{y}+p(x,y)u+q(x,y)\quad (1)$

The natural approach, after our experience with ordinary differential equations and with simple algebraic equations, is attempt a factorisation. Let's see how for this takes us.

We would expect factoring the left hand of (1) to give us an equivalent equation of the form

$a(x,y)(D_{x}+\alpha _{+}(x,y)D_{y})(D_{x}+\alpha _{-}(x,y)D_{y})u\,$

and we can immediately divide through by a. This suggests that those particular combinations of first order derivatives will play a special role.

Now, when studying first order PDE's we saw that such combinations were equivalent to the derivatives along characteristic curves. Effectively, we changed to a coordinate system defined by the characteristic curve and the initial curve.

Here, we have two combinations of first order derivatives each of which may define a different characteristic curve. If so, the two sets of characteristics will define a natural coordinate system for the problem, much as in the first order case.

In the new coordinates we will have

$D_{x}+\alpha _{+}(x,y)D_{y}=D_{r}\quad D_{x}+\alpha _{-}(x,y)D_{y}=D_{s}$

with each of the factors having become a differentiation along its respective characteristic curve, and the left hand side will become simply u_rs giving us an equation of the form

$u_{rs}=A(r,s)u_{r}+B(r,s)u_{s}+C(r,s)u+D(r,s)\,$

If A, B, and C all happen to be zero, the solution is obvious. If not, we can hope that the simpler form of the left hand side will enable us to make progress.

However, before we can do all this, we must see if (1) can actually be factored.

Multiplying out the factors gives

$u_{xx}+{\frac {b(x,y)}{a(x,y)}}u_{xy}+{\frac {c(x,y)}{a(x,y)}}u_{yy}=u_{xx}+(\alpha _{+}+\alpha _{-})u_{xy}+\alpha _{+}\alpha _{-}u_{yy}$

On comparing coefficients, and solving for the α's we see that they are the roots of

$a(x,y)\alpha ^{2}+b(x,y)\alpha +c(x,y)=0\,$

Since we are discussing real functions, we are only interested in real roots, so the existence of the desired factorization will depend on the discriminant of this quadratic equation.

• ${\mbox{If }}b(x,y)^{2}>4a(x,y)c(x,y)$

then we have two factors, and can follow the procedure outlined above. Equations like this are called hyperbolic

• ${\mbox{If }}b(x,y)^{2}=4a(x,y)c(x,y)$

then we have only factor, giving us a single characteristic curve. It will be natural to use distance along these curves as one coordinate, but the second must be determined by other considerations.

The same line of argument as before shows that use the characteristic curve this way gives a second order term of the form u_rr, where we've only taken the second derivative with respect to one of the two coordinates. Equations like this are called parabolic

• ${\mbox{If }}b(x,y)^{2}<4a(x,y)c(x,y)$

then we have no real factors. In this case the best we can do is reduce the second order terms to the simplest possible form satisfying this inequality, i.e u_rr+u_ss

It can be shown that this reduction is always possible. Equations like this are called elliptic

It can be shown that, just as for first order PDEs, discontinuities propagate along characteristics. Since elliptic equations have no real characteristics, this implies that any discontinuities they may have will be restricted to isolated points; i.e., that the solution is almost everywhere smooth.

This is not true for hyperbolic equations. Their behavior is largely controlled by the shape of their characteristic curves.

These differences mean different methods are required to study the three types of second equation. Fortunately, changing variables as indicated by the factorisation above lets us reduce any second order PDE to one in which the coefficients of the second order terms are constant, which means it is sufficient to consider only three standard equations.

$u_{xx}+u_{yy}=0\quad u_{xx}-u_{yy}=0\quad u_{xx}-u_{y}=0$

We could also consider the cases where the right hand side of these equations is a given function, or proportional to u or to one of its first order derivatives, but all the essential properties of hyperbolic, parabolic, and elliptic equations are demonstrated by these three standard forms.

While we've only demonstrated the reduction in two dimensions, a similar reduction applies in higher dimensions, leading to a similar classification. We get, as the reduced form of the second order terms,

$a_{1}{\frac {\partial ^{2}u}{\partial x_{1}^{2}}}+a_{2}{\frac {\partial ^{2}u}{\partial x_{2}^{2}}}+\cdots +a_{n}{\frac {\partial ^{2}u}{\partial x_{n}^{2}}}$

where each of the a_is is equal to either 0, +1, or -1.

If all the a_is have the same sign the equation is elliptic

If any of the a_is are zero the equation is parabolic

If exactly one of the a_is has the opposite sign to the rest the equation is hyperbolic

In 2 or 3 dimensions these are the only possibilities, but in 4 or more dimensions there is a fourth possibility: at least two of the a_is are positive, and at least two of the a_is are negative.

Such equations are called ultrahyperbolic. They are less commonly encountered than the other three types, so will not be studied here.

When the coefficients are not constant, an equation can be hyperbolic in some regions of the xy plane, and elliptic in others. If so, different methods must be used for the solutions in the two regions.

Elliptic

Standard form, Laplace's equation:

$\nabla ^{2}h=0$

Quote equation in spherical and cylindrical coordinates, and give full solution for cartesian and cylindrical coordinates. Note averaging property Comment on physical significance, rotation invariance of laplacian.

Hyperbolic

Standard form, wave equation:

$\nabla ^{2}h=c^{2}h_{tt}$

Solution, any sum of functions of the form

$h=f(\mathbf {k} \cdot \mathbf {x} -\omega t)\quad \omega =|\mathbf {k} |c$

These are waves. Compare with solution from separating variables. Domain of dependence, etc.

Parabolic

The canonical parabolic equation is the diffusion equation:

$\nabla ^{2}h=h_{t}$

Here, we will consider some simple solutions of the one-dimensional case.

The properties of this equation are in many respects intermediate between those of hyperbolic and elliptic equation.

As with hyperbolic equations but not elliptic, the solution is well behaved if the value is given on the initial surface t=0.

However, the characteristic surfaces of this equation are the surfaces of constant t, thus there is no way for discontinuities to propagate to positive t.

Therefore, as with elliptic equations but not hyberbolic, the solutions are typically smooth, even when the initial conditions aren't.

Furthermore, at a local maximum of h, its Laplacian is negative, so h is decreasing with t, while at local minima, where the Laplacian will be positive, h will increase with t. Thus, initial variations in h will be smoothed out as t increases.

In one dimension, we can learn more by integrating both sides,

${\begin{matrix}\int _{-a}^{b}h_{t}dt&=&\int _{-a}^{b}h_{xx}dx\\{\frac {d}{dt}}\int _{-a}^{b}h\,dt&=&\left[h_{x}\right]_{-a}^{b}\end{matrix}}$

Provided that h_x tends to zero for large x, we can take the limit as a and b tend to infinity, deducing

${\frac {d}{dt}}\int _{-\infty }^{\infty }h\,dt$

so the integral of h over all space is constant.

This means this PDE can be thought of as describing some conserved quantity, initially concentrated but spreading out, or diffusing, over time.

This last result can be extended to two or more dimensions, using the theorems of vector calculus.

We can also differentiate any solution with respect to any coordinate to obtain another solution. E.g. if h is a solution then

$\nabla ^{2}h_{x}=\partial _{x}\nabla ^{2}h=\partial _{x}\partial _{t}h=\partial _{t}h_{x}$

so h_x also satisfies the diffusion equation.

Similarity solution

Looking at this equation, we might notice that if we make the change of variables

$r=\alpha x\quad \tau =\alpha ^{2}t$

then the equation retains the same form. This suggests that the combination of variables x²/t, which is unaffected by this variable change, may be significant.

We therefore assume this equation to have a solution of the special form

$h(x,t)=f(\eta ){\mbox{ where }}\eta ={\frac {x}{t^{1 \over 2}}}$

then

$h_{x}=\eta _{x}f_{\eta }=t^{-{1 \over 2}}f_{\eta }\quad h_{t}=\eta _{t}f_{\eta }=-{\frac {\eta }{2t}}f_{\eta }$

and substituting into the diffusion equation eventually gives

$f_{\eta \eta }+{\frac {\eta }{2}}f_{\eta }=0$

which is an ordinary differential equation.

Integrating once gives

$f_{\eta }=Ae^{-{\frac {\eta ^{2}}{4}}}$

Reverting to h, we find

${\begin{matrix}h_{x}&=&{\frac {A}{\sqrt {t}}}e^{-{\frac {\eta ^{2}}{4}}}\\h&=&{\frac {A}{\sqrt {t}}}\int _{-\infty }^{x}e^{-s^{2}/4t}ds+B\\&=&A\int _{-\infty }^{x/2{\sqrt {t}}}e^{-z^{2}}dz+B\end{matrix}}$

This last integral can not be written in terms of elementary functions, but its values are well known.

In particular the limiting values of h at infinity are

$h(-\infty ,t)=B\quad h(\infty ,t)=B+A{\sqrt {\pi }},$

taking the limit as t tends to zero gives

$h=\left\{{\begin{matrix}B&x<0\\B+A{\sqrt {\pi }}&x>0\end{matrix}}\right.$

We see that the initial discontinuity is immediately smoothed out. The solution at later times retains the same shape, but is more stretched out.

The derivative of this solution with respect to x

$h_{x}={\frac {A}{\sqrt {t}}}e^{-x^{2}/4t}$

is itself a solution, with h spreading out from its initial peak, and plays a significant role in the further analysis of this equation.

The same similarity method can also be applied to some non-linear equations.

Separating variables

We can also obtain some solutions of this equation by separating variables.

$h(x,t)=X(x)T(t)\Rightarrow X^{\prime \prime }T=X{\dot {T}}$

giving us the two ordinary differential equations

${\frac {d^{2}X}{dx^{2}}}+k^{2}X=0\quad {\frac {dT}{dt}}=-kT$

and solutions of the general form

$h(x,t)=Ae^{-kt}\sin(kx+\alpha )\,$