Advanced Mathematics for Engineers and Scientists/Printable Version

Advanced Mathematics for Engineers and Scientists

The current, editable version of this book is available in Wikibooks, the open-content textbooks collection, at
https://en.wikibooks.org/wiki/Advanced_Mathematics_for_Engineers_and_Scientists

Permission is granted to copy, distribute, and/or modify this document under the terms of the Creative Commons Attribution-ShareAlike 3.0 License.

The Front Cover

Nomenclature

ODE	Ordinary Differential Equation
PDE	Partial Differential Equation
BC	Boundary Condition
IVP	Initial Value Problem
BVP	Boundary Value Problem
IBVP	Initial Boundary Value Problem

Common Operators

Operators are shown applied to the scalar $u(x_{1},x_{2},\cdots ,x_{n})$ or the vector field $\mathbf {v} (x_{1},x_{2},\cdots ,x_{n})=(v_{1},v_{2},\cdots ,v_{n})\,$.

Notation	Common names and other notation	Description and notes	Definition in Cartesian coordinates
${\frac {\partial u}{\partial x_{i}}}$	Partial derivative, $u_{x_{i}},\ \partial _{x_{i}}u\,$	The rate of change of $u$ with respect to $x_{i}$, holding the other independent variables constant.	$\lim _{\Delta x_{i}\to 0}{\frac {u(x_{1},\cdots ,x_{i}+\Delta x_{i},\cdots ,x_{n})-u}{\Delta x_{i}}}$
${\frac {du}{dx_{i}}}$	Derivative, total derivative, ${\frac {\mathrm {d} u}{\mathrm {d} x_{i}}}\,$	The rate of change of $u$ with respect to $x_{i}$. If $u$ is multivariate, this derivative will typically depend on the other variables following a path.	${\frac {\partial u}{\partial x_{1}}}{\frac {dx_{1}}{dx_{i}}}+\cdots +{\frac {\partial u}{\partial x_{n}}}{\frac {dx_{n}}{dx_{i}}}$
$\nabla u$	Gradient, del operator, $\mathrm {grad} \ u\,$	Vector that describes the direction and magnitude of the greatest rate of change of a function of more than one variable. The symbol $\nabla$ is called nabla.	$\left({\frac {\partial u}{\partial x_{1}}},\cdots ,{\frac {\partial u}{\partial x_{n}}}\right)$
$\nabla ^{2}u$	Laplacian, Scalar Laplacian, Laplace operator, $\Delta u,\ (\nabla \cdot \nabla )u\,$	A measure of the concavity of $u$, equivalently a comparison of the value of $u$ at some point to neighboring values.	${\frac {\partial ^{2}u}{\partial x_{1}^{2}}}+\cdots +{\frac {\partial ^{2}u}{\partial x_{n}^{2}}}$
$\nabla \cdot \mathbf {v}$	Divergence, $\mathrm {div} \ \mathbf {v} \,$	A measure of "generation", in other words how much the vector field acts as a source or sink at a point.	${\frac {\partial v_{1}}{\partial x_{1}}}+\cdots +{\frac {\partial v_{n}}{\partial x_{n}}}$
$\nabla \times \mathbf {v}$	Curl, rotor, circulation density, $\mathrm {curl} \ \mathbf {v} ,\ \mathrm {rot} \ \mathbf {v} \,$	A vector that describes the rate of rotation of a (normally 3D) vector field and the corresponding axis of rotation.	$\left({\frac {\partial v_{3}}{\partial x_{2}}}-{\frac {\partial v_{2}}{\partial x_{3}}},{\frac {\partial v_{1}}{\partial x_{3}}}-{\frac {\partial v_{3}}{\partial x_{1}}},{\frac {\partial v_{2}}{\partial x_{1}}}-{\frac {\partial v_{1}}{\partial x_{2}}}\right)$
$\nabla ^{2}\mathbf {v}$	Vector Laplacian	Similar to the (scalar) Laplacian. Note however, that it is generally not equal to the element-by-element Laplacian of a vector.	$\nabla (\nabla \cdot \mathbf {\mathbf {v} } )-\nabla \times (\nabla \times \mathbf {\mathbf {v} } )$

3D Operators in Different Coordinate Systems

Cartesian representations appear in the table above. The $(r,\theta ,\phi )=(\mathrm {distance,azimuth,colatitude} )$ convention is used for spherical coordinates.

Operator	Cylindrical	Spherical
$\nabla u$	$\left({\frac {\partial u}{\partial r}},{\frac {1}{r}}{\frac {\partial u}{\partial \theta }},{\frac {\partial u}{\partial z}}\right)\,$	$\left({\frac {\partial u}{\partial r}},{\frac {1}{r\sin(\phi )}}{\frac {\partial u}{\partial \theta }},{\frac {1}{r}}{\frac {\partial u}{\partial \phi }}\right)\,$
$\nabla ^{2}u$	${\frac {1}{r}}{\frac {\partial }{\partial r}}\left(r{\frac {\partial u}{\partial r}}\right)+{\frac {1}{r^{2}}}{\frac {\partial ^{2}u}{\partial \theta ^{2}}}+{\frac {\partial ^{2}u}{\partial z^{2}}}\,$	${\frac {1}{r^{2}}}{\frac {\partial }{\partial r}}\left(r^{2}{\frac {\partial u}{\partial r}}\right)+{\frac {1}{r^{2}\sin(\phi )}}{\frac {\partial ^{2}u}{\partial \theta ^{2}}}+{\frac {1}{r^{2}\sin(\phi )}}{\frac {\partial }{\partial \phi }}\left(\sin(\phi ){\frac {\partial u}{\partial \phi }}\right)\,$
$\nabla \cdot \mathbf {v}$	${\frac {1}{r}}{\frac {\partial }{\partial r}}\left(rv_{r}\right)+{\frac {1}{r}}{\frac {\partial v_{\theta }}{\partial \theta }}+{\frac {\partial v_{z}}{\partial z}}\,$	${\frac {1}{r^{2}}}{\frac {\partial }{\partial r}}\left(r^{2}v_{r}\right)+{\frac {1}{r\sin(\phi )}}{\frac {\partial v_{\theta }}{\partial \theta }}+{\frac {1}{r\sin(\phi )}}{\frac {\partial }{\partial \phi }}\left(\sin(\phi )v_{\phi }\right)\,$

Introduction to Partial Differential Equations

This book is intended as a Partial Differential Equations (PDEs) reference for individuals who already possess a firm understanding of ordinary differential equations and have at least a basic idea of what a partial derivative is.

This book is meant to be easily readable to engineers and scientists while still being (almost) interesting enough for mathematics students. Be advised that in-depth proofs of such matters as series convergence, uniqueness, and existence will not be given; this fact will appall some and elate others. This book is meant more toward solving or at the very least extracting information out of problems involving partial differential equations. The first few chapters are built to be especially simple to understand so that, say, the interested engineering undergraduate can benefit; however, later on important and more mathematical topics such as vector spaces will be introduced and used.

What follows is a quick intro for the uninitiated, with analogies to ordinary differential equations.

What is a Partial Differential Equation?

Ordinary differential equations (ODEs) arise naturally whenever a rate of change of some entity is known. This may be the rate of increase of a population, the rate of change of velocity, or maybe even the rate at which soldiers die on a battlefield. ODEs describe such changes of discrete entities. Respectively, this may be the capita of a population, the velocity of a particle, or the size of a military force.

More than one entity may be described with more than one ODE. For example, cloth is very often simulated in computer graphics as a grid of particles interconnected by springs, with Newton's law (an ODE) applied to each "cloth particle". In three dimensions, this would result in 3 second order ODEs written and solved for each particle.

Partial differential equations (PDEs) are analogous to ODEs in that they involve rates of change; however, they differ in that they treat continuous media. For example, the cloth could just as well be considered to be some kind of continuous sheet. This approach would most likely lead to only 3 (maybe 4) partial differential equations, which would represent the entire continuous sheet, instead of a set of ODEs for each particle.

This continuum approach is a very different way of looking at things. It may or may not be favorable: in the case of cloth, the resulting PDE system would be too difficult to solve, and so the computer graphics industry goes with a particle based approach (but a prime counterexample is a fluid, which would be represented by a PDE system most of the time).
While PDEs may not be straightforward to solve on a computer, they have a major advantage over ODEs when applicable: it is nearly impossible to gain any analytical insight from a huge system of particles, while a relatively small PDE system can reveal much insight, even if it won't yield an analytic solution.

But PDEs don't strictly describe continuum mechanics. As with anything mathematical, they are what you make of them.

The Character of Partial Differential Equations

The solution of an ODE can be represented as a function of one variable. For example, the position of the Earth may be represented by coordinates with respect to, say, the sun, and each of these coordinates would be functions of time. Note that the effects of other celestial bodies would certainly affect the solution, but it would still be expressible strictly as a function of time.

The solution of a PDE will, in general, depend on more than one variable. An example is a vibrating string: the deflection of the string will depend both on time and which part of the string you're looking at.

The solution of an ODE is called a trajectory. It may be represented graphically by one or more curves. The solution of a PDE, however, could be a surface, a volume, or something else, depending on how many variables are involved and how they're interpreted.

In general, PDEs are complicated to solve. Concepts such as separation of variables or integral transformations tend to work very differently. One significant difficulty is that the solution of a PDE depends very strongly on the initial/boundary conditions (ICs/BCs). An ODE typically yields a general solution, which involves one or more constants which may be determined from one or more ICs/BCs. PDEs, however, do not easily yield such general solutions. A solution method that works for one initial boundary value problem (IBVP) may be useless for a different IBVP.

PDEs tend to be more difficult to solve numerically as well. Most of the time, an ODE can be expressed in terms of its highest order derivative, and can be solved on a computer very easily with knowledge of the ICs (boundary value problems are a little more complicated), using well established and more or less generally applicable methods, such as Runge Kutta (RK). With this in mind, an ODE may be solved quickly by entering the equation and its ICs/BCs into the right application and pressing the solve button. An IBVP for a PDE, however, will typically require its own specialized solution, and it may take much effort to make the solution more than, say, second order accurate.

An Early Example

Many of the concepts of the previous section may be summarized in this example. We won't deal with the PDE just yet.

Consider heat flow along a laterally insulated rod. In other words, the heat is only flowing along the rod but not into the surrounding air. Let's call the temperature of the rod $u$, and let $u=u(x,t)$, where $t$ is time and $x$ represents the position along the rod. As the temperature depends both on time and position along the rod, this is exactly what $u=u(x,t)$ says. It is the change of the heat distribution over time. See the graphic below to get an idea.
Let's say that the rod has unitless length $1$, and that its initial temperature (again unitless) is known to be $u(x,0)=\sin(\pi x)$. This states the initial condition, which depends on $x$. The function is a simple hump between 0 and 1. Check for yourself with maxima (http://maxima.sourceforge.net or on android): plot2d(sin(x*%pi),[x,0,1])
Let's also say that the temperature is somehow fixed to $0$ at both ends of the rod, i.e. at $x=0$ and at $x=1$. This would result in $u(0,t)=u(1,t)=0$, which specifies boundary conditions. The BCs state that for all t, $u=0$ at $x=0$ and $x=1$.

A PDE can be written to describe the situation. This and the IC/BCs form an initial boundary value problem (IBVP). The solution to this IBVP is (with a physical constant taken to be $1$):

$u(x,t)=e^{-\pi ^{2}t}\sin(\pi x)\,$

Note that:

$u(x,0)=e^{-\pi ^{2}0}\sin(\pi x)=\sin(\pi x)\qquad {\mbox{(it satisfies the IC)}}\,$

$u(0,t)=e^{-\pi ^{2}t}\sin(\pi \cdot 0)=0\,$

$u(1,t)=e^{-\pi ^{2}t}\sin(\pi \cdot 1)=0\qquad {\mbox{(it satisfies the BCs)}}\,$

It also satisfies the PDE, but (again) that'll come later.

This solution may be interpreted as a surface, it's shown in the figure below with $x$ going from $0$ to $1$, and $t$ going from $0$ to $0.5$. That is, the distribution of heat is changing over time as the heat flows and dissipates.

$u(x,t)$ from $t=0$ to $t=0.5$ and $x=0$ to $x=1$.

Surfaces may or may not be the best way to convey information, and in this case a possibly better way to draw the picture would be to graph $u(x,t)$ as a curve at several different choices of $t$, this is portrayed below.

$u(x,t)$ in the domain of interest for various interesting values of $t$.

PDEs are extremely diverse, and their ICs and BCs can radically affect their solution method. As a result, the best (read: easiest) way to learn is by looking at many different problems and how they're solved.

Parallel Plate Flow: Easy IC

Parallel Plate Flow: Easy IC

Formulation

As with ODEs, separation of variables is easy to understand and works well whenever it works. For ODEs, we use the substitution rule to allow antidifferentiation, but for PDEs it's a very different process involving letting dependencies pass through the partial derivatives.

A fluid mechanics example will be used.

Consider two plates parallel to each other of huge extent, separated by a distance of 1. Fluid is smoothly flowing between these two plates in only one direction (call it x). This may be seen in the picture below.

Visualization of the parallel plate flow problem.

After some assumptions, the following PDE may be obtained to describe the fluid flow:

${\frac {\partial u}{\partial t}}=\nu {\frac {\partial ^{2}u}{\partial y^{2}}}-{\frac {P_{x}}{\rho }}\,$

This linear PDE is the result of simplifying the Navier-Stokes equations, a large nonlinear PDE system which describes fluid flow. u is the velocity of the fluid in the x direction, ρ is the density of the fluid, ν is the kinematic viscosity (metaphorically speaking, how hard the molecules of the fluid rub against each other), and P_x is the pressure gradient or pressure gradient vector field. Note that u = u(y, t), there is no dependence on x. In other words, the state of the fluid upstream is no different from the state downstream. Notice that we do not consider turbulences and that the state of the flow varies between the upper plate and the lower plate as the speed u is 0 close to the plates.

u(y, t) is a velocity profile. Fluid mechanics typically is concerned with velocity fields, contrary to rigid body mechanics in which the position of an object is what is important. In other words, with a rigid object all 'points' of the object move at the same speed. With a fluid we get a velocity field and each point is moving at its own speed.

The ratio P_x/ρ describes the driving force; it's a pressure change (gradient) along the x direction. If P_x is negative, then the pressure downstream (positive x) is smaller than the pressure upstream (negative x) and the fluid will flow left to right, i.e., u(y, t) will generally be positive.

Now on to create a specific problem: let's say that a constant negative pressure gradient was applied for a long time, until the velocity profile was steady (steady means "not changing with time"). Then the pressure gradient is suddenly removed, and without this driving force the fluid will slow down and stop. That is the assumption we are going to use for our example calculations. We assume the flow to be steady and then decaying as the pressure (expressed as P_x/ρ) is removed. Hence we can remove the term -P_x/ρ from our model as well, see (PDE) below.

Initial flow profile.

Let's say that before the pressure was removed, the velocity profile was u(y, t) = sin(π y). With velocity profile we mean, as the velocities are measured on a cross section of the flow, they form a hump. This would make sense: the friction dictates less motion near the plates (see next paragraph), so we could expect a maximum velocity near the centerline (y = 1/2). This assumed profile isn't really correct, but will serve as an example for now. It's graphed at right in the domain of interest.

Before getting into the math, one more thing is needed: boundary conditions. In this case, the BC is called the no slip condition, which states that the velocity of a fluid at a wall (boundary) is equal to the velocity of the wall. If we weren't making this assumption, the fluid would be moving like a rigid object. Since the velocities of the walls (or plates) in this problem are both zero, the velocity of the fluid must be zero at these two boundaries. The BCs are then u(0, t) = 0 (bottom plate) and u(1, t) = 0 (top plate).

The IBVP is:

${\frac {\partial u}{\partial t}}=\nu {\frac {\partial ^{2}u}{\partial y^{2}}}\qquad {\mbox{(PDE)}}\,$

$u(y,0)=\sin(\pi y)\qquad {\mbox{(IC)}}\,$

$u(0,t)=0\,$

$u(1,t)=0\qquad {\mbox{(BCs)}}\,$

Notice that the initial condition IC is the same hump from the first section of this book, just turned sideways. This IC implies that the flow is already flowing when we start our calculations. We are kind of calculating the decay of the flow speed.

Separation

Variables are separated the following way: we assume that $u(y,t)=Y(y)T(t)$, where Y and T are (unknown) functions respectively of y and t. This form is substituted into the PDE:

${\frac {\partial }{\partial t}}(u(y,t))=\nu {\frac {\partial ^{2}}{\partial y^{2}}}(u(y,t))\,$

Using $u(y,t)=Y(y)T(t)$ yields:

${\frac {\partial }{\partial t}}(Y(y)T(t))=\nu {\frac {\partial ^{2}}{\partial y^{2}}}(Y(y)T(t))\,$

$Y(y){\frac {\partial }{\partial t}}(T(t))=\nu T(t){\frac {\partial ^{2}}{\partial y^{2}}}(Y(y))\,$

$Y(y)T'(t)=\nu T(t)Y''(y)\,$

${\frac {T'(t)}{\nu T(t)}}={\frac {Y''(y)}{Y(y)}}\,$

Look carefully at the last equation: the left side of the equation depends strictly on t, and the right side strictly on y, and they are equal. t may be varied independently of y and they'd still be equal, and y may be varied independently of t and they'd still be equal. This can only happen if both sides are constant. This may be shown as follows:

${\frac {d}{dt}}\left({\frac {T'(t)}{\nu T(t)}}\right)={\frac {d}{dt}}\left({\frac {Y''(y)}{Y(y)}}\right)\,$

${\frac {d}{dt}}\left({\frac {T'(t)}{\nu T(t)}}\right)=0\,$

Taking the derivative of both sides, turns the right-hand-side into a constant, 0. The left-hand-side is left as is. Then taking the integrals of both sides yields:

${\frac {T'(t)}{\nu T(t)}}={\mbox{a constant}}\,$

Integration of the ordinary derivative recovers the left side but leaves the right side a constant. It follows by similarity that Y''/Y is a constant as well.

The constant in question is called the separation constant. We could simply give it a letter, such as A, but a good choice of the constant will make work easier later. In this case the best choice is -k². This will be justified later (but it should be reemphasized that it may be notated any way you want, assuming it can span the domain).

${\frac {T'(t)}{\nu T(t)}}={\frac {Y''(y)}{Y(y)}}=-k^{2}\,$

${\Big \Downarrow }$

${\frac {T'(t)}{\nu T(t)}}=-k^{2}\,$

${\frac {Y''(y)}{Y(y)}}=-k^{2}\,$

The variables are now separated. The last two equations are two ODEs which may be solved independently (in fact, the Y equation is an eigenvalue problem), though they both contain an unknown constant. Note that ν was kept for the T equation. This choice makes the solution slightly easier, but is again completely arbitrary.

${\frac {T'(t)}{\nu T(t)}}=-k^{2}\,$

Rearrange and note that the notion buried in the expression below is that a function is equal to its own derivative, which kind of strikes the Euler number bell.

$T'(t)=-k^{2}\nu T(t)\,$

Then solve.

$T(t)=C_{1}e^{-k^{2}\nu t}\,$

Here is a solution plucked from Ted Woollett's 'Maxima by Example', Chapter 3, '3.2.3 Exact Solution Using desolve'. Maxima's output is not shown but should be easily reconcilable with what has been written.

(%i1) de:(-%k^2*v*T(t))-('diff(T(t),t)) = 0;

(%i2) gsoln:desolve(de,T(t));

The same goes for our right-hand-side.

${\frac {Y''(y)}{Y(y)}}=-k^{2}\,$

$Y''(y)=-k^{2}Y(y)\,$

And solve.

$Y(y)=C_{2}\cos(ky)+C_{3}\sin(ky)\,$

In Maxima the solution goes like this. When Maxima is asking for 'zero' or 'nonzero' then type 'nonzero;':

(%i1)de:(-%k^2*Y(y))-('diff(Y(y),y,2));

(%i2)atvalue('diff(Y(y),y), y=0, Y(0)*%k )$

(%i3)gsoln:desolve(de,Y(y));

The overall solution will be $u(y,t)$, still with unknown constants, and as of now the product of Y and T. So we are plugging the partial solutions for $Y(y)$ and $T(t)$ back in:

$u(y,t)=Y(y)T(t)=C_{1}e^{-k^{2}\nu t}(C_{2}\cos(ky)+C_{3}\sin(ky))$

$u(y,t)=e^{-k^{2}\nu t}(A\cos(ky)+B\sin(ky))$

Note that C₁ has been multiplied into C₂ and C₃, reducing the number of arbitrary constants. Because an unknown constant multiplied by an unknown constant yields still an unknown constant.

The IC or BCs should now be applied. If the IC was applied first, coefficients would be equated and all of the constants would be determined. However, the BCs may or may not have been fulfilled (in this case they would, but you're not generally so lucky). So to be safe, the BCs will be applied first:

$u(0,t)=e^{-k^{2}\nu t}(A\cos(k\cdot 0)+B\sin(k\cdot 0))=0$

$A\cdot 1+B\cdot 0=0\Rightarrow A=0\qquad {\mbox{(no slip along lower boundary)}}\,$

So A being zero eliminates the $\cos(k\cdot 0)$ term.

$u(1,t)=e^{-k^{2}\nu t}B\sin(k\cdot 1)=0\,$

$B\sin(k\cdot 1)=0\Rightarrow B=0\quad {\mbox{or}}\quad k=n\pi \qquad {\mbox{(no slip along upper boundary)}}\,$

If we took B = 0, the solution would have just been u(y, t) = 0 (often called the trivial solution), which would satisfy the BCs and the PDE but couldn't possibly satisfy the IC. So, we take k = nπ instead, where n is any integer. After applying the BCs we have:

Decaying flow.

$u(y,t)=e^{-(n\pi )^{2}\nu t}B\sin(n\pi y)\,$

Then we need to apply the IC to it. Per the IC from above $u(y,0)$ is:

$u(y,0)=\sin(\pi y)\,$

Per the BC $u(y,0){\mbox{with t = 0}}$ is, see above. Setting those two definitions of $u(y,0)$ equal is:

$e^{0}B\sin(n\pi y)=\sin(\pi y)\,$

Since t = 0 as per the IC assumption, $e^{-(n\pi )^{2}\nu t}$ is becoming $e^{0}$. The equality can only hold if B = 1 and n = 1, allowing us to simply remove those two constants from the function in its BC incarnation. That's it! The complete solution is:

$u(y,t)=e^{-\pi ^{2}\nu t}\sin(\pi y)\,$

It's worth verifying that the IC, BCs, and PDE are all satisfied by this. Also notice that the solution is a product of a function of t and a function of y. The graph at the right is illustrating this. Observe that the profile is plotted for different values of νt, rather than specifying some ν and graphing different values for t. Remember that v is the kinematic viscosity. Hence the decay of flow speed also depends on it. Looking at the solution, t and ν appear only once and they're multiplying, so it's natural to do this. A dimensionless time could have been introduced from the beginning.

So what happens? The fluid starts with its initial profile and slows down exponentially. Note that with x replaced with y and t replaced with νt, this is exactly the same as the result of heat flow in a rod as shown in the introduction. This is not a coincidence: the PDE for the rod describes diffusion of heat, the PDE for the parallel plates describes diffusion of momentum.

Take a second look at the separation constant, -k². The square is convenient, without it the solution for Y(y) would have involved square roots. Without the negative sign, the solution would have involved exponentials instead of sinusoids, so the constant would have come out imaginary.

The assumption that u(y, t) = Y(y)T(t) is justified by the physics of the problem: it would make sense that the profile would keep its general shape (due to Y(y)), only it'd get flattened over time as the fluid slows down (due to T(t)).

Quickly recap what we did. First we took a model PDE from Navier-Stokes and simplified it by making an assumption about the pressure being removed. Then we sort of first-stage solved it by separating Y and T. Then we applied the BC (boundary conditions) and the IC (initial condition) yielding our final solution.

Parallel Plate Flow: Realistic IC

Parallel Plate Flow: Realistic IC

The Steady State

The initial velocity profile chosen in the last problem agreed with intuition but honestly came out of thin air. A more realistic development follows.

The problem stated that (to come up with an IC) the fluid was under a pressure difference for some time, so that the flow became steady aka flowing steadily. "Steady" is another way of saying "not changing with time", and "not changing with time" is another way of saying that:

${\frac {\partial u}{\partial t}}=0\,$

Putting this into the PDE from the previous section:

${\frac {\partial u}{\partial t}}=\nu {\frac {\partial ^{2}u}{\partial y^{2}}}-{\frac {P_{x}}{\rho }}\,$

${\Big \Downarrow }$

$0=\nu {\frac {\partial ^{2}u}{\partial y^{2}}}-{\frac {P_{x}}{\rho }}\quad \Rightarrow \quad {\frac {P_{x}}{\nu \rho }}={\frac {d^{2}u}{dy^{2}}}\,$

Independent of $t$, the PDE became an ODE with variables separated and thus we can integrate.

${\frac {d^{2}u}{dy^{2}}}={\frac {P_{x}}{\nu \rho }}\quad \Rightarrow \quad u={\frac {P_{x}}{2\nu \rho }}y^{2}+C_{1}y+C_{0}\,$

The no slip condition results in the following BCs: $u=0$ at $y=0$ and $y=1$. We can plug the BC values into the integrated ODE and resolve the Cs.

$0={\frac {P_{x}}{2\nu \rho }}\cdot 0^{2}+C_{1}\cdot 0+C_{0}\quad \Rightarrow \quad C_{0}=0\qquad {\mbox{(Bottom plate BC: u = 0, y = 0)}}\,$

$0={\frac {P_{x}}{2\nu \rho }}\cdot 1^{2}+C_{1}\cdot 1\quad \Rightarrow \quad C_{1}=-{\frac {P_{x}}{2\nu \rho }}\qquad {\mbox{(Top plate BC: u = 0, y = 1)}}\,$

Inserting the Cs and and simplifying yields:

$u={\frac {P_{x}}{2\nu \rho }}(y^{2}-y)\,$

For the sake of example, take $P_{x}/(2\nu \rho )=-4$ (recall that a negative pressure gradient causes left to right flow). Also note that this is a constant gradient or slope. This gives a parabola which starts at $0$, increases to a maximum of $1$ at $y=1/2$, and returns to $0$ at $y=1$.

This parabola looks pretty much identical to the sinusoid previously used (you must zoom in to see a difference). However, even more so on the narrow domain of interest, the two are very different functions (look at their taylor expansions, for example). Using the parabola instead of the sine function results in a much more involved solution.

So this derives the steady state flow, which we will use as an improved, realistic IC. Recall that the problem is about a fluid that's initially in motion that is coming to a stop due to the absence of a driving force. The IBVP (Initial Boundary Value Problem) is now subtly different:

${\frac {\partial u}{\partial t}}=\nu {\frac {\partial ^{2}u}{\partial y^{2}}}\qquad {\mbox{(PDE)}}\,$

$u(y,0)=4y-4y^{2}\qquad {\mbox{(IC)}}\,$

$u(0,t)=0\,$

$u(1,t)=0\qquad {\mbox{(BCs)}}\,$

Separation

Since the only difference from the problem in the last section is the IC, the variables may be separated and the BCs applied with no difference, giving:

$u(y,t)=e^{-(n\pi )^{2}\nu t}B\sin(n\pi y)\,$

But now we're stuck (after applying the BCs)! Applying the IC makes the $e^{-(n\pi )^{2}\nu t}$ term go away as t = 0, which is the IC. However, then the IC function can't be made to match:

$u(y,0)=4y-4y^{2}\,$

$B\sin(n\pi y)=4y-4y^{2}\,$

What went wrong? It was the assumption that $u(y,t)=Y(y)T(t)$. The fact that the IC couldn't be fulfilled means that the assumption was wrong. It should be apparent now why the IC was chosen to be $\sin(\pi y)$ in the previous section.

We can proceed however, thanks to the linearity of the problem. Another detour is necessary, it gets long.

Linearity (the superposition principle specifically) says that if $u_{1}$ is a solution to the BVP (not the whole IBVP, only the BVP, Boundary Value Problem, the BCs applied) and so is another $u_{2}$, then a linear combination, $C_{1}u_{1}+C_{2}u_{2}$, is also a solution.

Let's take a step back and suppose that the IC was

$u(y,0)=\sin(\pi y)+1/5\sin(3\pi y).$

This is no longer a realistic flow problem but it contains the first two terms of what is called a Fourier sine expansion, see these examples of Fourier sine expansions. We are going to generalize this below. Let's now use this expression and equate it to the half way solution (BCs applied) with $e^{-(n\pi )^{2}\nu t}$ being eliminated as t = 0:

$B\sin(n\pi y)=\sin(\pi y)+{\frac {1}{5}}\sin(3\pi y)\,$

And it still can't match. However, observe that the individual terms in the IC can. We simply set the constants to values making both sides match:

$B_{1}\sin(n_{1}\pi y)=\sin(\pi y)\Rightarrow n_{1}=1\ {\mbox{and}}\ B_{1}=1\,$

$B_{3}\sin(n_{3}\pi y)={\frac {1}{5}}\sin(3\pi y)\Rightarrow n_{3}=3\ {\mbox{and}}\ B_{3}={\frac {1}{5}}\,$

Note the subscripts are used to identify each term: they reflect the integer $n$ from the separation constant. Solutions may be obtained for each individual term of the IC, identified with $n$:

$u_{1}(y,t)=e^{-(1\pi )^{2}\nu t}1\sin(1\pi y)=e^{-\pi ^{2}\nu t}\sin(\pi y)$        $n=1\ {\mbox{and}}\ B=1\,$

$u_{3}(y,t)=e^{-(3\pi )^{2}\nu t}{\frac {1}{5}}\sin(3\pi y)=e^{-9\pi ^{2}\nu t}{\frac {1}{5}}\sin(3\pi y)\,$        $n=3\ {\mbox{and}}\ B={\frac {1}{5}}\,$

Linearity states that the sum of these two solutions is also a solution to the BVP (no need for new constants):

$u_{1+3}(y,t)=e^{-\pi ^{2}\nu t}\sin(\pi y)+{\frac {1}{5}}e^{-9\pi ^{2}\nu t}\sin(3\pi y)\,$

So we added the solutions and got a new solution... what is this good for? Try setting $t=0$:

$u_{1+3}(y,0)=e^{-\pi ^{2}\nu \cdot 0}\sin(\pi y)+{\frac {1}{5}}e^{-9\pi ^{2}\nu \cdot 0}\sin(3\pi y)=\sin(\pi y)+{\frac {1}{5}}\sin(3\pi y)\,$

Each component solution satisfies the BVP, and the sum of these just happened to satisfy our surrogate IC. The IBVP with IC $u(y,0)=\sin(\pi y)+1/5\sin(3\pi y)$ is now solved. It would work the same way for any linear combination of sine functions whose half frequencies are $n\pi$. "Linear combination" means a sum of terms, each multiplied by a constant. The sum is assumed to converge and be term by term differentiable.

Let's do what we just did in a more generalized fashion. First, we make our IC a linear combination of sines (with $e^{-(n\pi )^{2}\nu t}$ eliminated as t = 0), in fact, infinitely many of them. But each successive term has to 'converge', it can't stray wildly all over the place.

$u(y,0)=\sum _{n=1}^{\infty }B_{n}\sin(n\pi y)\qquad {\mbox{(IC: an arbitrary linear combination of sines)}}\,$

Second, find the n and B for each term assuming t = 0 (the IC), then plug them back into each term making no assumptions about t, leaving t as is.

$u_{n}(y,t)=e^{-(n\pi )^{2}\nu t}B_{n}\sin(n\pi y)\qquad {\mbox{(solution of}}\ n^{th}\ {\mbox{term)}}\,$

Third, sum up all the terms with their individual n and Bs.

$u(y,t)=\sum _{n=1}^{\infty }u_{n}(y,t)=\sum _{n=1}^{\infty }e^{-(n\pi )^{2}\nu t}B_{n}\sin(n\pi y)\qquad {\mbox{(sum of solutions)}}\,$

Fourth, plug t = 0 into the sum of terms and recover the IC from the first step.

$u(y,0)=\sum _{n=1}^{\infty }e^{-(n\pi )^{2}\nu \cdot 0}B_{n}\sin(n\pi y)=\sum _{n=1}^{\infty }B_{n}\sin(n\pi y)\qquad {\mbox{(IC recovered)}}\,$

So we went full circle on this example but found the n and Bs because we were able to equate/satisfy each term with the IC. Now we can solve the problem if the IC is a linear combination of sine functions. But the IC for this problem isn't such a sum, it's just a stupid parabola. Or is it?

Series Construction

In the 19th century, a man named Joseph Fourier took a break from helping Napoleon take over the world to ask an important question while studying this same BVP (concerning heat flow): can a function be expressed as a sum of sinusoids, similar to a taylor series? The short answer is yes, if a few reasonable conditions apply as we have already indicated. The long answer follows, and this section is a longer answer.

A function meeting certain criteria may indeed be expanded into a sum of sines, cosines, or both. In our case, all that is needed to accomplish this expansion is to find the coefficients $B_{n}$. A little trick involving an integral makes this possible.

The sine function has a very important property called orthogonality. There are many flavors of this, which will be served in the next chapter. Relevant to this problem is the following:

$\int _{0}^{1}2\sin(m\pi y)\sin(n\pi y)\,dy={\begin{cases}1,&m=n\\0,&m\neq n\end{cases}}$

A quick hint may help. Orthogonality literally means two lines at a right angle to each other. These lines could be vectors, each with its own tuple of coordinates. If those two vectors are at a right angle to each other, multiplying and summing their coordinate tuples always yields zero (in Euclidean space). The method of multiplying and summing is also used to determine whether two functions are orthogonal. Using this definition, our multiplied and integrated functions above are orthogonal most of the time, but not always.

Let's call the IC $\phi (y)$ to generalize it. We equate the IC with its expansion, meaning the linear combination of sines, and then apply some craftiness. And remember that our goal is to reproduce a parabolic function from linearly combined sines:

$\sum _{n=1}^{\infty }B_{n}\sin(n\pi y)=\phi (y)\,$

$2\sin(m\pi y)\cdot \sum _{n=1}^{\infty }B_{n}\sin(n\pi y)=2\sin(m\pi y)\cdot \phi (y)\,$

$\sum _{n=1}^{\infty }B_{n}\cdot 2\sin(m\pi y)\sin(n\pi y)=2\sin(m\pi y)\phi (y)\,$

$\int _{0}^{1}\sum _{n=1}^{\infty }B_{n}\cdot 2\sin(m\pi y)\sin(n\pi y)dy=\int _{0}^{1}2\sin(m\pi y)\phi (y)\ dy\,$

$\sum _{n=1}^{\infty }B_{n}\int _{0}^{1}2\sin(m\pi y)\sin(n\pi y)dy=\int _{0}^{1}2\sin(m\pi y)\phi (y)\ dy\,$

$B_{m}=\int _{0}^{1}2\sin(m\pi y)\phi (y)\ dy\,$

In the last step, all of the terms in the sum became $0$ except for the $m^{th}$ term where $m=n$, the only case where we get $1$ for the otherwise orthogonal sine functions. This isolates and explicitly defines $B_{m}$ which is the same as $B_{n}$ as m = n. The expansion for $\phi (y)$ is then:

$\phi (y)=\sum _{m=1}^{\infty }(\int _{0}^{1}2\sin(m\pi y)\phi (y)\ dy)\cdot \sin(m\pi y)\,$

Or equivalently:

$\phi (y)=\sum _{m=1}^{\infty }B_{m}\sin(m\pi y)\qquad ;\ B_{m}=\int _{0}^{1}2\sin(m\pi y)\phi (y)\ dy\,$

Many important details have been left out for later in a devoted chapter; one noteworthy detail is that this expansion is only approximating the parabola (very superficially) on the interval $0\leq y\leq 1$, not say from $-\infty$ to $+\infty$.

This expansion may finally be combined with the sum of sines solution to the BVP developed previously. Note that the last equation looks very similar to $u(y,0)$. Following from this:

$u(y,0)=\sum _{n=1}^{\infty }B_{n}\sin(n\pi y)\,$

$u(y,0)=\sum _{n=1}^{\infty }\int _{0}^{1}2\sin(n\pi y)\phi (y)\ dy\cdot \sin(n\pi y)\qquad \,$

So the expansion will satisfy the IC given as $\phi (y)$ (surprised?). The full solution for the problem with arbitrary IC is then:

$u(y,t)=\sum _{n=1}^{\infty }e^{-(n\pi )^{2}\nu t}B_{n}\sin(n\pi y)\,$

$u(y,t)=\sum _{n=1}^{\infty }e^{-(n\pi )^{2}\nu t}\int _{0}^{1}2\sin(n\pi y)\phi (y)\ dy\cdot \sin(n\pi y)\,$

In this problem specifically, the IC is $\phi (y)=4y-4y^{2}$, so:

$B_{n}=\int _{0}^{1}2\sin(n\pi y)(4y-4y^{2})dy={\frac {8}{n^{3}\pi ^{3}}}(2-2\cos(n\pi )-n\pi \sin(n\pi ))\,$

Sines and cosines appear from the integration dependent only on $n\pi$. Since $n$ is an integer, these can be made more aesthetic.

$\sin(n\pi )=0\qquad ;\ n{\mbox{ is an integer.}}\,$

$\cos(n\pi )=(-1)^{n}\qquad ;\ n{\mbox{ is an integer.}}\,$

$\qquad {\Big \Downarrow }$

$B_{n}={\frac {16-16(-1)^{n}}{n^{3}\pi ^{3}}}\,$

Note that for even $n$, $B_{n}=0$. Putting everything together finally completes the solution to the IBVP:

$u(y,t)=\sum _{n=1}^{\infty }e^{-(n\pi )^{2}\nu t}{\frac {16-16(-1)^{n}}{n^{3}\pi ^{3}}}\sin(n\pi y)\,$

There are many interesting things to observe. To begin with, $u(y,t)$ is not a product of a function of $y$ and a function of $t$. Such a solution was assumed in the beginning, proved to be wrong, but eventually happened to yield a solution anyway thanks to linearity and what is called a Fourier sine expansion.

A careful look at the procedure reveals something that may be disturbing: this lengthy solution is strictly valid for the given BCs. Thanks to the definition of $\phi (y)$, the solution is generic as far as the IC is concerned (the IC doesn't even need to match the BCs), however a slight change in either BC would mandate starting over almost from the beginning.

The parabolic IC, which looks very similar to the sine function used in the previous section, is wholly to blame (or thank once you understand the beauty of a Fourier series!) for the infinite sum. It is interesting to approximate the first several numeric values of the sequence $B_{n}$:

${\begin{aligned}B_{1}&\approx 1.03205\\B_{3}&\approx 0.03822\\B_{5}&\approx 0.00826\\B_{7}&\approx 0.00301\,\end{aligned}}$

Recall that the even terms are all $0$. The first term by far dominates, this makes sense since the first term already looks very, very similar to the parabola. Recall that $n^{2}$ appears in an exponential, making the higher terms even smaller for time not too close to $0$.

Change of Variables

Change of Variables

As with ODEs (ordinary differential equations), a PDE (partial differential equation, or more accurately, the initial-boundary value problem (IBVP) as a whole) may be made more amenable with the help of some kind of modification of variables. So far, we've dealt only with boundary conditions that specify the value of u, which represented fluid velocity, as zero at the boundaries. Though fluid mechanics can get more complicated than that (understatement of the millennium), let's look at heat transfer now for the sake of variety.

As hinted previously, the one dimensional diffusion equation can also describe heat flow in one dimension. Think of how heat could flow in one dimension: one possibility is a rod that's completely laterally insulated, so that the heat will flow only along the rod and not across it (be aware, though, it is possible to consider heat loss/gain along the rod without going two dimensional).

If this rod has finite length, heat could flow in and out of the uninsulated ends. A 1D rod can have at most two ends (it can also have one or zero: the rod could be modeled as "very long"), and the boundary conditions could specify what happens at these ends. For example, the temperature could be specified at a boundary, or maybe the flow of heat, or maybe some combination of the two.

The equation for heat flow is usually given as:

${\frac {\partial u}{\partial t}}=\alpha {\frac {\partial ^{2}u}{\partial x^{2}}}\,$

Which is the same as the equation for parallel plate flow, only with ν replaced with α and y replaced with x.

Fixed Temperatures at Boundaries

Let's consider a rod of length 1, with temperatures specified (fixed) at the boundaries. The IBVP is:

${\frac {\partial u}{\partial t}}=\alpha {\frac {\partial ^{2}u}{\partial x^{2}}}\,$

$u(x,0)=\varphi (x)\,$

$u(0,t)=u_{0}\,$

$u(1,t)=u_{1}\,$

φ(x) is the temperature at t = 0. Look at what the BCs say: For all time, the temperature at x = 0 is u₀ and at x = 1 is u₁. Note that this could be just as well a parallel plate problem: u₀ and u₁ would represent wall velocities.

The PDE is easily separable, in basically the same way as in previous chapters:

${\frac {\partial }{\partial t}}(u(x,t))=\alpha {\frac {\partial ^{2}}{\partial x^{2}}}(u(x,t))\,$

${\frac {\partial }{\partial t}}(X(x)T(t))=\alpha {\frac {\partial ^{2}}{\partial x^{2}}}(X(x)T(t))\,$

$X(x){\frac {\partial }{\partial t}}(T(t))=\alpha T(t){\frac {\partial ^{2}}{\partial x^{2}}}(X(x))\,$

$X(x)T'(t)=\alpha T(t)X''(x)\,$

${\frac {T'(t)}{\alpha T(t)}}={\frac {X''(x)}{X(x)}}\,$

${\Big \Downarrow }$

${\frac {T'(t)}{\alpha T(t)}}={\frac {X''(x)}{X(x)}}=-k^{2}\,$

${\Big \Downarrow }$

${\frac {T'(t)}{\alpha T(t)}}=-k^{2}\Rightarrow T'(t)=-k^{2}\alpha T(t)\Rightarrow T(t)=C_{1}e^{-k^{2}\alpha t}\,$

${\frac {X''(x)}{X(x)}}=-k^{2}\Rightarrow X''(x)=-k^{2}X(x)\Rightarrow X(x)=C_{2}\cos(kx)+C_{3}\sin(kx)\,$

${\Big \Downarrow }$

$u(x,t)=X(x)T(t)=C_{1}e^{-k^{2}\alpha t}(C_{2}\cos(kx)+C_{3}\sin(kx))$

$u(x,t)=e^{-k^{2}\alpha t}(A\cos(kx)+B\sin(kx))$

Now, substitute the BCs:

$u_{0}=u(0,t)\,$

$u_{0}=e^{-k^{2}\alpha t}(A\cos(k\cdot 0)+B\sin(k\cdot 0))\,$

$u_{0}=Ae^{-k^{2}\alpha t}\,$

$u_{1}=u(1,t)\,$

$u_{1}=e^{-k^{2}\alpha t}(A\cos(k\cdot 1)+B\sin(k\cdot 1))\,$

$u_{1}=e^{-k^{2}\alpha t}(A\cos(k)+B\sin(k))\,$

We can't proceed. Among other things, the presence of t in the exponential factor (previously divided out) prevents anything from coming out of this.

This is another example of the fact that the assumption that u(x, t) = X(x)T(t) was wrong. The only thing that prevents us from getting a solution would be the non-zero BCs. This is where changing variables will help: a new variable v(x, t) will be defined in terms of u which will be separable.

Think of how v(x, t) could be defined to make its BCs zero ("homogeneous"). One way would be:

$u(x,t)=v(x,t)+h(x)\,$

This form is inspired from the appearance of the BCs, and it can be readily seen:

$u(0,t)=u_{0}\,$

$u(1,t)=u_{1}\,$

${\Big \Downarrow }$

$v(0,t)+h(0)=u_{0}\,$

$v(1,t)+h(1)=u_{1}\,$

If h(0) = u₀ and h(1) = u₁, v(x, t) would indeed have zero BCs. Pretty much any choice of h(x) satisfying these conditions would do it, but only one is the best choice. Making the substitution into the PDE:

${\frac {\partial }{\partial t}}(u(x,t))=\alpha {\frac {\partial ^{2}}{\partial x^{2}}}(u(x,t))\,$

${\frac {\partial }{\partial t}}(v(x,t)+h(x))=\alpha {\frac {\partial ^{2}}{\partial x^{2}}}(v(x,t)+h(x))\,$

${\frac {\partial v}{\partial t}}=\alpha {\frac {\partial ^{2}v}{\partial x^{2}}}+\alpha {\frac {\partial ^{2}h}{\partial x^{2}}}\,$

So now the PDE has been messed up by the new term involving h. This will thwart separation...

...unless that last term happens to be zero. Rather then hoping it's zero, we can demand it (the best choice hinted above), and put the other requirements on h(x) next to that:

${\frac {d^{2}h}{dx^{2}}}=0\,$

$h(0)=u_{0}\,$

$h(1)=u_{1}\,$

Note that the partial derivative became an ordinary derivative since h is a function of x only. The above constitutes a pretty simple boundary value problem, with unique solution:

$h(x)=(u_{1}-u_{0})x+u_{0}\,$

It's just a straight line. Note that this is what would arise if the steady state (time independent) problem were solved for u(x). In other words, h could've been pulled out of one's ass readily just looking at the physics of the situation.

The problem now reduces to finding v(x, t). The IBVP for this would be:

${\frac {\partial v}{\partial t}}=\alpha {\frac {\partial ^{2}v}{\partial x^{2}}}\,$

$v(x,0)=\varphi (x)-h(x)=\varphi (x)-((u_{1}-u_{0})x+u_{0})\,$

$v(0,t)=0\,$

$v(1,t)=0\,$

Note that the IC changed under the transformation. The solution to this IBVP was found in a past chapter through separation of variables and superposition to be:

$v(x,t)=\sum _{n=1}^{\infty }e^{-(n\pi )^{2}\alpha t}\int _{0}^{1}2\sin(n\pi x)(\varphi (x)-((u_{1}-u_{0})x+u_{0}))\ dx\cdot \sin(n\pi x)\,$

u(x, t) may now be found simply by adding h(x), according to how the variable change was defined:

$u(x,t)=v(x,t)+h(x)\,$

$u(x,t)=\sum _{n=1}^{\infty }e^{-(n\pi )^{2}\alpha t}\int _{0}^{1}2\sin(n\pi x)(\varphi (x)-h(x))\ dx\cdot \sin(n\pi x)+(u_{1}-u_{0})x+u_{0}\,$

This solution looks like the sum of a steady state portion (that's h(x)) and a transient portion (that's v(x)):

Visualization of the change of variables.

Time Varying Temperatures at Boundaries

Note that this wouldn't work so nicely with non-constant BCs. For example, if the IBVP were:

${\frac {\partial u}{\partial t}}=\alpha {\frac {\partial ^{2}u}{\partial x^{2}}}\,$

$u(x,0)=\varphi (x)\,$

$u(0,t)=u_{0}(t)\,$

$u(1,t)=u_{1}(t)\,$

Then, transforming it would require h = h(x, t). Reusing u(x, t) = v(x, t) + h(x, t) introduced previously would eventually lead to:

${\frac {\partial v}{\partial t}}+{\frac {\partial h}{\partial t}}=\alpha {\frac {\partial ^{2}v}{\partial x^{2}}}+\alpha {\frac {\partial ^{2}h}{\partial x^{2}}}\,$

$v(x,0)=\varphi (x)-h(x,0)\,$

$v(0,t)=0\,$

$v(1,t)=0\,$

Where, to simplify the PDE above:

${\frac {\partial h}{\partial t}}=\alpha {\frac {\partial ^{2}h}{\partial x^{2}}}\,$

$h(x,0)={\mbox{anything}}\,$

$h(0,t)=u_{0}(t)\,$

$h(1,t)=u_{1}(t)\,$

Which doesn't really make anything simpler, despite freedom in the choice of IC.

But this isn't completely useless. Note that the PDE for h was chosen to simplify the PDE for v(x, t) (would lead to the terms involving h to cancel out), which may lead to the question: Was this necessary?

The answer is no. If that were the case, the PDE we picked for h would not be satisfied, and that would result in extra terms in the PDE for v(x, t). The no-longer-separable IBVP for v(x, t) could, however, be solved via an eigenfunction expansion, whose full story will be told later. It's worth noting though, that an eigenfunction expansion would require homogenous BCs, so the transformation was necessary.

So this problem has to be put aside without any conclusion for now. I told you that BCs can mess everything up.

Pressure Driven Transient Parallel Plate Flow

Now back to fluid mechanics. Previously, we dealt with flow that was initially moving but slowing down because of resistance and the absence of a driving force. Maybe, it'd be more interesting if we had a fluid initially at rest (ie, zero IC) but set into motion by some constant pressure difference. The IBVP for such a case would be:

${\frac {\partial u}{\partial t}}=\nu {\frac {\partial ^{2}u}{\partial y^{2}}}-{\frac {P_{x}}{\rho }}\,$

$u(y,0)=0\,$

$u(0,t)=0\,$

$u(1,t)=0\,$

This PDE with the pressure term was described previously. That pressure term is what drives the flow; it is assumed constant.

The intent of the change of variables would be to remove the pressure term from the PDE (which prevents separation) while keeping the BCs homogeneous.

One path to take would be to add something to u(x, t), either a function of t or a function of y, so that differentiation would leave behind a constant that could cancel the pressure term out. Adding a function of t would be very unfavorable since it'd result in time dependent BCs, so let's try a function of y:

$u(y,t)=v(y,t)+f(y)\,$

Substituting this into the PDE:

${\frac {\partial }{\partial t}}(v(y,t)+f(y))=\nu {\frac {\partial ^{2}}{\partial y^{2}}}(v(y,t)+f(y))-{\frac {P_{x}}{\rho }}\,$

${\frac {\partial v}{\partial t}}=\nu {\frac {\partial ^{2}v}{\partial y^{2}}}+\nu {\frac {\partial ^{2}f}{\partial y^{2}}}-{\frac {P_{x}}{\rho }}\,$

This procedure will simplify the PDE and preserve the BCs only if the following conditions hold:

$0=\nu {\frac {d^{2}f}{dy^{2}}}-{\frac {P_{x}}{\rho }}\,$

$f(0)=0\,$

$f(1)=0\,$

The first condition, an ODE, is required to simplify the PDE for v(y, t), it will result in cancellation of the last two terms. The other two conditions are chosen to preserve the homogeneous BCs of the problem (note that if the BCs of u(y, t) weren't homogeneous, the BCs on f(y) would need to be picked to amend that).

The solution to the BVP above is simply:

$f(y)={\frac {P_{x}}{2\rho \nu }}(y^{2}-y)\,$

So f(y) was successfully determined. Note that this function is symmetric about y = 1/2. The IBVP for v(y, t) becomes:

${\frac {\partial v}{\partial t}}=\nu {\frac {\partial ^{2}v}{\partial y^{2}}}\,$

$v(y,0)=-f(y)={\frac {P_{x}}{2\rho \nu }}(y-y^{2})\,$

$v(0,t)=0\,$

$v(1,t)=0\,$

This is the same IBVP we've been beating to death for some time now. The solution for v(y, t) is:

$v(y,t)=\sum _{n=1}^{\infty }e^{-(n\pi )^{2}\alpha t}\int _{0}^{1}2\sin(n\pi y)\left({\frac {P_{x}}{2\rho \nu }}(y-y^{2})\right)\ dy\cdot \sin(n\pi y)\,$

And the solution for u(y, t) follows from how the variable change was defined:

$u(y,t)=v(y,t)+f(y)\,$

$u(y,t)=\sum _{n=1}^{\infty }e^{-(n\pi )^{2}\nu t}\int _{0}^{1}2\sin(n\pi y)\left({\frac {P_{x}}{2\rho \nu }}(y-y^{2})\right)\ dy\cdot \sin(n\pi y)+{\frac {P_{x}}{2\rho \nu }}(y^{2}-y)\,$

$u(y,t)={\frac {P_{x}}{2\rho \nu }}(y^{2}-y)-{\frac {P_{x}}{2\rho \nu }}\sum _{n=1}^{\infty }e^{-(n\pi )^{2}\nu t}\cdot {\frac {4(-1)^{n}-4}{n^{3}\pi ^{3}}}\sin(n\pi y)\,$

This solution fits what we expect: it starts flat and approaches the parabolic profile quickly. This is the same parabola derived as the steady state flow in the realistic IC chapter; the integral was evaluated for integer n, simplifying it.

A careful look at the solution reveals something interesting: this is just decaying parallel plate flow "in reverse". Instead of the flow starting parabolic and gradually approaching u = 0, it starts with u = 0 and gradually approaches a parabola.

Time Dependent Diffusivity

In this example we'll change time, an independent variable, instead of changing the dependent variable. Consider the following IBVP:

${\frac {\partial u}{\partial t}}={\frac {1}{1+t^{2}}}{\frac {\partial ^{2}u}{\partial x^{2}}}\,$

$u(x,0)=\varphi (x)\,$

$u(0,t)=0\,$

$u(1,t)=0\,$

Note that this is a separable; a transformation isn't really necessary, however it'll be easier since we can reuse past solutions if it can be transformed into something familiar.

Let's not get involved with the physics of this and just call it a diffusion problem. It could be diffusion of momentum (as in fluid mechanics), diffusion of heat (heat transfer), diffusion of a chemical (chemistry), or simply a mathematician's toy. In other words, a confession: it was purposely made up to serve as an example.

The (time dependent) factor in front of the second derivative is called the diffusivity. Previously, it was a constant α (called "thermal diffusivity") or constant ν ("kinematic viscosity"). Now, it decays with time.

To simplify the PDE via a transformation, we look for ways in which the factor could cancel out. One way would be to define a new time variable, call it τ and leave it's relation to t arbitrary. The chain rule yields:

${\frac {\partial u}{\partial t}}={\frac {\partial u}{\partial \tau }}\cdot {\frac {d\tau }{dt}}\,$

Substituting this into the PDE:

${\frac {\partial u}{\partial \tau }}\cdot {\frac {d\tau }{dt}}={\frac {1}{1+t^{2}}}{\frac {\partial ^{2}u}{\partial x^{2}}}\,$

Note now that the variable t will completely disappear (divide out in this case) from the equation if:

${\frac {d\tau }{dt}}={\frac {1}{1+t^{2}}}\Rightarrow \tau =\arctan(t)+C\,$

C is completely arbitrary. However, the best choice of C is the one that makes τ = 0 when t = 0, since this wouldn't change the IC which is defined at t = 0; so, take C = 0. Note that the BCs wouldn't change either way, unless they were time dependent, in which they would change no matter what C is chosen. The IBVP is turned into:

${\frac {\partial u}{\partial \tau }}={\frac {\partial ^{2}u}{\partial x^{2}}}\,$

$u(x,0)=\varphi (x)\,$

$u(0,\tau )=0\,$

$u(1,\tau )=0\,$

Digging up the solution and restoring the original variable:

$u(x,\tau )=\sum _{n=1}^{\infty }e^{-(n\pi )^{2}\tau }\int _{0}^{1}2\sin(n\pi x)\varphi (x)\ dx\cdot \sin(n\pi x)\,$

${\Big \Downarrow }$

$u(x,t)=\sum _{n=1}^{\infty }e^{-(n\pi )^{2}\arctan(t)}\int _{0}^{1}2\sin(n\pi x)\varphi (x)\ dx\cdot \sin(n\pi x)\,$

Note that, unlike any of the previous examples, the physics of the problem (if there were any) couldn't have helped us. It's also worth mentioning that the solution doesn't limit to u = 0 for long time.

Concluding Remarks

Changing variables works a little differently for PDEs in the sense that you have a lot of freedom thanks to partial differentiation. In this chapter, we picked what seemed to be a good general form for the transformation (inspired by whatever prevented easy solution), wrote down a bunch of requirements, and defined the transformation to uniquely satisfy the requirements. Doing the same for ODEs can often degrade to a monkey with typewriter situation.

Many simple little changes go without saying. For example, we've so far worked with rods of length "1" or plates separated by a distance of "1". What if the rod was 5 m long? Then space would have to be nondimensionalized using the following transformation:

$x=5\ {\mbox{m}}\cdot {\hat {x}}\,$

Simple nondimensionalization is, well, simple; however for PDEs with more terms it can lead to scale analysis which can lead to perturbation theory which will all have to be explained in a later chapter.

It's worth noting that the physics of the IBVP very often suggest what kind of transformation needs to be done. Even some nonlinear problems can be solved this way.

This topic isn't nearly over, changes of variables will be dealt with again in future chapters.

The Laplacian and Laplace's Equation

The Laplacian and Laplace's Equation

By now, you've most likely grown sick of the one dimensional transient diffusion PDE we've been playing with:

${\frac {\partial u}{\partial t}}=\alpha {\frac {\partial ^{2}u}{\partial x^{2}}}\,$

Make no mistake: we're not nearly done with this stupid thing; but for the sake of variety let's introduce a fresh new equation and, even though it's not strictly a separation of variables concept, a really cool quantity called the Laplacian. You'll like this chapter; it has many pretty pictures in it.

Graph of $u(x)=x^{3}-x$.

The Laplacian

The Laplacian is a linear operator in Euclidean n-space. There are other spaces with properties different from Euclidean space. Note also that operator here has a very specific meaning. As a function is sort of an operator on real numbers, our operator is an operator on functions, not on the real numbers. See here for a longer explanation.

We'll start with the 3D Cartesian "version". Let $u=u(x,y,z,t)$. The Laplacian of the function $u$ is defined and notated as:

$\nabla ^{2}u=\Delta u=\operatorname {div} \ \operatorname {grad} \ u={\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}+{\frac {\partial ^{2}u}{\partial z^{2}}}\,$

So the operator is taking the sum of the nonmixed second derivatives of $u$ with respect to the Cartesian space variables $x,y,z$. The "del squared" notation is preferred since the capital delta can be confused with increments and differences, and ${\mbox{div grad }}u$ is too long and doesn't involve pretty math symbols. The Laplacian is also known as the Laplace operator or Laplace's operator, not to be confused with the Laplace transform. Also, note that if we had only taken the first partial derivatives of the function $u$, and put them into a vector, that would have been the gradient of the function $u$. The Laplacian takes the second unmixed derivatives and adds them up.

In one dimension, recall that the second derivative measures concavity. Suppose $y=f(x)$; if $f''(x)$ is positive, $y$ is concave up, and if $f''(x)$ is negative, $y$ is concave down, see the graph below with the straight up or down arrows at various points of the curve. The Laplacian may be thought of as a generalization of the concavity concept to multivariate functions.

This idea is demonstrated at the right, in one dimension: $u(x)=x^{3}-x$. To the left of $x=0$, the Laplacian (simply the second derivative here) is negative, and the graph is concave down. At $x=0$, the curve inflects and the Laplacian is $0$. To the right of $x=0$, the Laplacian is positive and the graph is concave up.

Concavity may or may not do it for you. Thankfully, there's another very important view of the Laplacian, with deep implications for any equation it shows itself in: the Laplacian compares the value of $u$ at some point in space to the average of the values of $u$ in the neighborhood of the same point. The three cases are:

• If $u$ is greater at some point than the average of its neighbors, $\scriptstyle \nabla ^{2}u<0$.
• If $u$ is at some point equal to the average of its neighbors, $\scriptstyle \nabla ^{2}u=0$.
• If $u$ is smaller at some point than the average of its neighbors, $\scriptstyle \nabla ^{2}u>0$.

So the laplacian may be thought of as, at some point $(x_{0},y_{0},z_{0})$:

$\nabla ^{2}u{\Big |}_{x=x_{0},y=y_{0},Z=Z_{0}}\approx {\mbox{average of u in the neighborhood of }}(x_{0},y_{0},z_{0})\ -\ u(x_{0},y_{0},z_{0},t)\,$

The neighborhood of $(x_{0},y_{0},z_{0})$.

The neighborhood of some point is defined as the open set that lies within some Euclidean distance δ (delta) from the point. Referring to the picture at right (a 3D example), the neighborhood of the point $(x_{0},y_{0},z_{0})$ is the shaded region which satisfies:

$(x-x_{0})^{2}+(y-y_{0})^{2}+(z-z_{0})^{2}<\delta ^{2}\,$

Note that our one dimensional transient diffusion equation, our parallel plate flow, involves the Laplacian:

${\frac {\partial u}{\partial t}}=\alpha {\frac {\partial ^{2}u}{\partial x^{2}}}=\alpha \nabla ^{2}u\,$

With this mentality, let's examine the behavior of this very important PDE. On the left is the time derivative and on the right is the Laplacian. This equation is saying that:

The rate of change of $u$ at some point is proportional to the difference between the average value of $u$ around that point and the value of $u$ at that point.

For example, if there's at some position a "hot spot" where $u$ is on average greater than its neighbors, the Laplacian will be negative and thus the time derivative will be negative, this will cause $u$ to decrease at that position, "cooling" it down. This is illustrated below. The arrows reflect upon the magnitude of the Laplacian and, by grace of the time derivative, the direction the curve will move.

Visualization of transient diffusion.

It's worth noting that in 3D, this equation fully describes the flow of heat in a homogeneous solid that's not generating it's own heat (like too much electricity through a narrow wire would).

Laplace's Equation

Laplace's equation describes a steady state condition, and this is what it looks like:

$\nabla ^{2}u=0\,$

Solutions of this equation are called harmonic functions. Some things to note:

• Time is absent. This equation describes a steady state condition.
• The absence of time implies the absence of an IC, so we'll be dealing with BVPs rather then IBVPs.
• In one dimension, this is the ODE of a straight line passing through the boundaries at their specified values.
• All functions that satisfy this equation in some domain are analytic (informally, an analytic function is equal to its Taylor expansion) in that domain.
• Despite appearances, solutions of Laplace's equation are generally not minimal surfaces.
• Laplace's equation is linear.

Laplace's equation is separable in the Cartesian (and almost any other) coordinate system. So, we shouldn't have too much problem solving it if the BCs involved aren't too convoluted.

Laplace's Equation on a Square: Cartesian Coordinates

Steady state conditions on a square.

Imagine a 1 x 1 square plate that's insulated top and bottom and has constant temperatures applied at its uninsulated edges, visualized to the right. Heat is flowing in and out of this thing steadily through the edges only, and since it's "thin" and "insulated", the temperature may be given as $u(x,y)$. This is the first time we venture into two spatial coordinates, note the absence of time.

Let's make up a BVP, referring to the picture:

${\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}=0\,$

${\begin{aligned}u(0,y)=0\qquad {\mbox{(Edge D)}}\\u(1,y)=0\qquad {\mbox{(Edge B)}}\\u(x,0)=0\qquad {\mbox{(Edge A)}}\\u(x,1)=1\qquad {\mbox{(Edge C)}}\end{aligned}}$

So we have one nonhomogeneous BC. Assume that $u(x,y)=X(x)Y(y)$:

${\frac {\partial ^{2}}{\partial x^{2}}}(X(x)Y(y))+{\frac {\partial ^{2}}{\partial y^{2}}}(X(x)Y(y))=0\,$

$X''(x)Y(y)+X(x)Y''(y)=0\,$

$-{\frac {X''(x)}{X(x)}}={\frac {Y''(y)}{Y(y)}}=k^{2}\,$

${\Big \Downarrow }$

$X''(x)=-k^{2}X(x)\,$

$Y''(y)=k^{2}Y(y)\,$

As with before, calling the separation constant $k^{2}$ in favor of just $k$ (or something) happens to make the problem easier to solve. Note that the negative sign was kept for the $X(x)$ equation: again, these choices happen to make things simpler. Solving each equation and combining them back into $u(x,y)$:

$X(x)=C_{1}\cos(kx)+C_{2}\sin(kx)\,$

$Y(y)=C_{3}e^{ky}+C_{4}e^{-ky}\,$

${\Big \Downarrow }$

$u(x,y)=X(x)Y(y)\,$

$u(x,y)=\left(C_{1}\cos(kx)+C_{2}\sin(kx)\right)\left(C_{3}e^{ky}+C_{4}e^{-ky}\right)\,$

At edge D:

$u(0,y)=0\,$

$\left(C_{1}\cos(k\cdot 0)+C_{2}\sin(k\cdot 0)\right)\left(C_{3}e^{ky}+C_{4}e^{-ky}\right)=0\,$

$C_{1}\left(C_{3}e^{ky}+C_{4}e^{-ky}\right)=0\Rightarrow C_{1}=0\,$

${\Big \Downarrow }$

$u(x,y)=C_{2}\sin(kx)\left(C_{3}e^{ky}+C_{4}e^{-ky}\right)\,$

Note that the constants can be merged, but we won't do it so that a point can be made in a moment. At edge A:

$u(x,0)=0\,$

$C_{2}\sin(kx)\left(C_{3}e^{k\cdot 0}+C_{4}e^{-k\cdot 0}\right)=0\,$

Taking $C_{2}$ as $0$ would satisfy this particular BC, however this would yield a plane solution of $u(x,y)=0$, which can't satisfy the temperature at edge C. This is why the constants weren't merged a few steps ago, to make it obvious that $C_{2}$ may not be $0$. So, we instead take $C_{3}=-C_{4}$ to satisfy the above, and then combine the three constants into one, call it $B$:

$u(x,y)=C_{2}\sin(kx)\left(C_{3}e^{ky}+C_{4}e^{-ky}\right)\qquad ;\ C_{3}=-C_{4}\,$

${\Big \Downarrow }$

$u(x,y)=B\sin(kx)\left(e^{ky}-e^{-ky}\right)\,$

Now look at edge B:

$u(1,y)=0\,$

$B\sin(k)\left(e^{ky}-e^{-ky}\right)=0\,$

It should go without saying by now that $B$ can't be zero, since this would yield $u(x,y)=0$ which couldn't satisfy the nonzero BC. Instead, we can take $k=n\pi$:

$u(x,y)=B\sin(n\pi x)\left(e^{n\pi y}-e^{-n\pi y}\right)\,$

As of now, this solution will satisfy 3 of the 4 BCs. All that is left is edge C, the nonhomogeneous BC.

$u(x,1)=1\,$

$B\sin(n\pi x)\left(e^{n\pi y}-e^{-n\pi y}\right)=1\,$

Neither $B$ nor $n$ can be contorted to fit this BC.

Since Laplace's equation is linear, a linear combination of solutions to the PDE is also a solution to the PDE. Another thing to note: since the BCs (so far) are homogeneous, we can add the solutions without worrying about nonzero boundaries adding up.

Though $u(x,y)$ as shown above will not solve this problem, we can try summing (based on $n$) solutions to form a linear combination which might solve the BVP as a whole:

$u(x,y)=\sum _{n=1}^{\infty }u_{n}(x,y)\,$

$u_{n}(x,y)=B_{n}\sin(n\pi x)\left(e^{n\pi y}-e^{-n\pi y}\right)\,$

${\Big \Downarrow }$

$u(x,y)=\sum _{n=1}^{\infty }B_{n}\sin(n\pi x)\left(e^{n\pi y}-e^{-n\pi y}\right)\,$

Assuming this form is correct (review Parallel Plate Flow: Realistic IC for motivation), let's again try applying the last BC:

$u(x,1)=1\,$

$\sum _{n=1}^{\infty }B_{n}\sin(n\pi x)\left(e^{n\pi \cdot 1}-e^{-n\pi \cdot 1}\right)=1\,$

It looks like it needs Fourier series methodology. Finding $B_{n}$ via orthogonality should solve this problem:

$2\sin(m\pi x)\cdot \sum _{n=1}^{\infty }B_{n}\sin(n\pi x)\left(e^{n\pi }-e^{-n\pi }\right)=2\sin(m\pi x)\cdot 1\,$

$\sum _{n=1}^{\infty }B_{n}\cdot 2\sin(m\pi x)\sin(n\pi x)\left(e^{n\pi }-e^{-n\pi }\right)=2\sin(m\pi x)\,$

$\int _{0}^{1}\sum _{n=1}^{\infty }B_{n}\cdot 2\sin(m\pi x)\sin(n\pi x)\left(e^{n\pi }-e^{-n\pi }\right)dx=\int _{0}^{1}2\sin(m\pi x)dx\,$

$\sum _{n=1}^{\infty }B_{n}\int _{0}^{1}2\sin(m\pi x)\sin(n\pi x)dx\left(e^{n\pi }-e^{-n\pi }\right)=2{\frac {1-\cos(m\pi )}{m\pi }}\,$

25 term partial sum of the series solution.

$B_{m}\left(e^{m\pi }-e^{-m\pi }\right)=2{\frac {1-\cos(m\pi )}{m\pi }}\,$

$B_{n}=2{\frac {1-(-1)^{n}}{n\pi \left(e^{n\pi }-e^{-n\pi }\right)}}\,$

$m$ was changed to $n$ in the last step. Also, for integer $n$, $\cos(n\pi )=(-1)^{n}$. Note that a Fourier sine expansion has been done. The solution to the BVP can finally be assembled:

$u(x,y)=\sum _{n=1}^{\infty }B_{n}\sin(n\pi x)\left(e^{n\pi y}-e^{-n\pi y}\right)\,$

$u(x,y)=\sum _{n=1}^{\infty }2{\frac {1-(-1)^{n}}{n\pi \left(e^{n\pi }-e^{-n\pi }\right)}}\sin(n\pi x)\left(e^{n\pi y}-e^{-n\pi y}\right)\,$

That solves it!

It's finally time to mention that the BCs are discontinuous at the points $(0,1)$ and $(1,1)$. As a result, the series should converge slowly at those points. This is clear from the plot at right: it's a 25 term partial sum (note that half of the terms are $0$), and it looks perfect except at $y=1$, especially near the discontinuities at $x=0$ and $x=1$.

Laplace's Equation on a Circle: Polar Coordinates

Now, we'll specify the value of $u$ on a circular boundary. A circle can be represented in Cartesian coordinates without too much trouble; however, it would result in nonlinear BCs which would render the approach useless. Instead, polar coordinates $(r,\theta )$ should be used, since in such a system the equation of a circle is very simple. In order for this to be realized, a polar representation of the Laplacian is necessary. Without going in to the details just yet, the Laplacian is given in (2D) polar coordinates:

$\nabla ^{2}u={\frac {\partial ^{2}u}{\partial r^{2}}}+{\frac {1}{r}}{\frac {\partial u}{\partial r}}+{\frac {1}{r^{2}}}{\frac {\partial ^{2}u}{\partial \theta ^{2}}}\,$

This result may be derived using differentials and the chain rule; it's not difficult but it's a little long. In these coordinates Laplace's equation reads:

${\frac {\partial ^{2}u}{\partial r^{2}}}+{\frac {1}{r}}{\frac {\partial u}{\partial r}}+{\frac {1}{r^{2}}}{\frac {\partial ^{2}u}{\partial \theta ^{2}}}=0\,$

Note that in going from Cartesian to polar coordinates, a price was paid: though still linear, Laplace's equation now has variable coefficients. This implies that after separation at least one of the ODEs will have variable coefficients as well.

Let's make up the following BVP, letting $u=u(r,\theta )$:

${\frac {\partial ^{2}u}{\partial r^{2}}}+{\frac {1}{r}}{\frac {\partial u}{\partial r}}+{\frac {1}{r^{2}}}{\frac {\partial ^{2}u}{\partial \theta ^{2}}}=0\,$

$u(1,\theta )=\phi (\theta )\,$

$u(r,\theta )\ {\mbox{is bounded inside the unit circle.}}\,$

This could represent a physical problem analogous to the previous one: replace the square plate with a disc. Note the apparent absence of sufficient BC to obtain a unique solution. The funny looking statement that u is bounded inside the domain of interest turns out to be the key to getting a unique solution, and it often shows itself in polar coordinates. It "makes up" for the "lack" of BCs. To separate, we as usual incorrectly assume that $u(r,\theta )=R(r)\Theta (\theta )$:

${\frac {\partial ^{2}}{\partial r^{2}}}(R(r)\Theta (\theta ))+{\frac {1}{r}}{\frac {\partial }{\partial r}}(R(r)\Theta (\theta ))+{\frac {1}{r^{2}}}{\frac {\partial ^{2}}{\partial \theta ^{2}}}(R(r)\Theta (\theta ))=0\,$

$R''(r)\Theta (\theta )+{\frac {1}{r}}R'(r)\Theta (\theta )+{\frac {1}{r^{2}}}R(r)\Theta ''(\theta )=0\,$

$r^{2}R''(r)+rR'(r)+R(r){\frac {\Theta ''(\theta )}{\Theta (\theta )}}=0\,$

$-{\frac {\Theta ''(\theta )}{\Theta (\theta )}}={\frac {r^{2}R''(r)+rR'(r)}{R(r)}}=k^{2}\,$

${\Big \Downarrow }$

$r^{2}R''(r)+rR'(r)-k^{2}R(r)=0\,$

$\Theta ''(\theta )+k^{2}\Theta (\theta )=0\,$

Once again, the way the negative sign and the separation constant are arranged makes the solution easier later on. These decisions are made mostly by trial and error.

The $R$ equation is probably one you've never seen before, it's a special case of the Euler differential equation (not to be confused with the Euler-Lagrange differential equation). There are a couple of ways to solve it, the most general method would be to change the variables so that an equation with constant coefficients is obtained. An easier way would be to note the pattern in the order of the coefficients and the order of the derivatives, and from there guess a power solution. Either way, the general solution to this simple case of Euler's ODE is given as:

$R(r)=C_{1}r^{k}+C_{2}r^{-k}\,$

This is a very good example problem since it goes to show that PDE problems very often turn into obscure ODE problems; we got lucky this time since the solution for $R$ was rather simple though its ODE looked pretty bad at first sight. The solution to the $\Theta (\theta )$ equation is:

$\Theta (\theta )=C_{3}\cos(k\theta )+C_{4}\sin(k\theta )\,$

Combining:

$u(r,\theta )=R(r)\Theta (\theta )\,$

$u(r,\theta )=\left(C_{1}r^{k}+C_{2}r^{-k}\right)(C_{3}\cos(k\theta )+C_{4}\sin(k\theta ))\,$

Now, this is where the English sentence condition stating that u must be bounded in the domain of interest may be invoked. As $\scriptstyle r\to 0$, the term involving $r^{-k}$ is unbounded. The only way to fix this is to take $C_{2}=0$. Note that if this problem were solved between two concentric circles, this term would be nonzero and very important. With that term gone, constants can be merged:

$u(r,\theta )=r^{k}(A\cos(k\theta )+B\sin(k\theta ))\,$

Only one condition remains: $u=\phi (\theta )$ on $r=1$, yet there are 3 constants. Let's say for now that:

$\phi (\theta )=2\cos(4\theta )+3\sin(4\theta )\,$

Then, it's a simple matter of equating coefficients to obtain:

$u(r,\theta )=r^{4}(2\cos(4\theta )+3\sin(4\theta ))\,$

Now, let's make the frequencies differ:

$\phi (\theta )=4\cos(3\theta )+\sin(10\theta )\,$

Equating coefficients won't work. However, if the IC were broken up into individual terms, the sum of the solution to the terms just happens to solve the BVP as a whole:

$\phi _{3}(\theta )=4\cos(3\theta )\quad ,\quad \phi _{10}(\theta )=\sin(10\theta )\,$

$u_{3}(r,\theta )=r^{3}4\cos(3\theta )\quad ,\quad u_{10}(r,\theta )=r^{10}\sin(10\theta )\,$

$u(r,\theta )=u_{3}(r,\theta )+u_{10}(r,\theta )=r^{3}4\cos(3\pi \theta )+r^{10}\sin(10\theta )\,$

Verify that the solution above is really equal to the BC at $r=1$:

$u(1,\theta )=1^{3}4\cos(3\theta )+1^{10}\sin(10\theta ))=4\cos(3\theta )+\sin(10\theta )=\phi (\theta )\,$

And, since Laplace's equation is linear, this must solve the PDE as well. What all of this implies is that, if some generic function $\phi (\theta )$ may be expressed as a sum of sinusoids with angular frequencies given by $n$, all that is needed is a linear combination of the appropriate sum. Notated:

$u(r,\theta )=\sum _{n=1}^{\infty }r^{n}(A_{n}\cos(n\theta )+B_{n}\sin(n\theta ))\,$

To identify the coefficients, substitute the BC:

$u(1,\theta )=\sum _{n=1}^{\infty }1^{n}(A_{n}\cos(n\theta )+B_{n}\sin(n\theta ))\,$

$\phi (\theta )=\sum _{n=1}^{\infty }(A_{n}\cos(n\theta )+B_{n}\sin(n\theta ))\,$

The coefficients $A_{n}$ and $B_{n}$ may be determined by a (full) Fourier expansion on $0\leq \theta \leq 2\pi$. Note that it's implied that $\phi (\theta )$ must have period $2\pi$ since we are solving this in a domain (a circle specifically) where $0\leq \theta \leq 2\pi$.

You probably don't like infinite series solutions. Well, it happens that through a variety of manipulations it's possible to express the full solution of this particular problem as:

$u(r,\theta )={\frac {1}{2\pi }}\int _{0}^{2\pi }\phi (\psi ){\frac {1-r^{2}}{1-2r\cos(\theta -\psi )+r^{2}}}\ d\psi \,$

This is called Poisson's integral formula.

Derivation of the Laplacian in Polar Coordinates

Though not necessarily a PDEs concept, it is very important for anyone studying this kind of math to be comfortable with going from one coordinate system to the next. What follows is a long derivation of the Laplacian in 2D polar coordinates using the multivariable chain rule and the concept of differentials. Know, however, that there are really many ways to do this.

Three definitions are all we need to begin:

$\nabla ^{2}u={\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}\,$

$x=r\cos(\theta )\,$

$y=r\sin(\theta )\,$

If it's known that $u=u(r,\theta )=u(r(x,y),\theta (x,y))$, then the chain rule may be used to express derivatives in terms of $r$ and $\theta$ alone. Two applications will be necessary to obtain the second derivatives. Manipulating operators as if they meant something on their own:

${\frac {\partial }{\partial x}}={\frac {\partial }{\partial r}}{\frac {\partial r}{\partial x}}+{\frac {\partial }{\partial \theta }}{\frac {\partial \theta }{\partial x}}\,$

Applying this to itself, treating the underlined bit as a unit dependent on $r$ and $\theta$:

${\frac {\partial ^{2}}{\partial x^{2}}}={\frac {\partial }{\partial x}}\left({\frac {\partial }{\partial x}}\right)={\frac {\partial }{\partial x}}\left({\underline {{\frac {\partial }{\partial r}}{\frac {\partial r}{\partial x}}+{\frac {\partial }{\partial \theta }}{\frac {\partial \theta }{\partial x}}}}\right)\,$

${\frac {\partial ^{2}}{\partial x^{2}}}={\frac {\partial }{\partial r}}\left({\underline {{\frac {\partial }{\partial r}}{\frac {\partial r}{\partial x}}+{\frac {\partial }{\partial \theta }}{\frac {\partial \theta }{\partial x}}}}\right){\frac {\partial r}{\partial x}}+{\frac {\partial }{\partial \theta }}\left({\underline {{\frac {\partial }{\partial r}}{\frac {\partial r}{\partial x}}+{\frac {\partial }{\partial \theta }}{\frac {\partial \theta }{\partial x}}}}\right){\frac {\partial \theta }{\partial x}}\,$

${\frac {\partial ^{2}}{\partial x^{2}}}=\left({\frac {\partial ^{2}}{\partial r^{2}}}{\frac {\partial r}{\partial x}}+{\frac {\partial }{\partial r}}{\frac {\partial ^{2}r}{\partial r\partial x}}+{\frac {\partial ^{2}}{\partial r\partial \theta }}{\frac {\partial \theta }{\partial x}}+{\frac {\partial }{\partial \theta }}{\frac {\partial ^{2}\theta }{\partial r\partial x}}\right){\frac {\partial r}{\partial x}}+\left({\frac {\partial ^{2}}{\partial \theta \partial r}}{\frac {\partial r}{\partial x}}+{\frac {\partial }{\partial r}}{\frac {\partial ^{2}r}{\partial \theta \partial x}}+{\frac {\partial ^{2}}{\partial \theta ^{2}}}{\frac {\partial \theta }{\partial x}}+{\frac {\partial }{\partial \theta }}{\frac {\partial ^{2}\theta }{\partial \theta \partial x}}\right){\frac {\partial \theta }{\partial x}}\,$

The above mess may be quickly simplified a little by manipulating the funny looking derivatives:

${\frac {\partial ^{2}r}{\partial r\partial x}}={\frac {\partial ^{2}}{\partial x\partial r}}(r)={\frac {\partial }{\partial x}}\left({\frac {\partial }{\partial r}}(r)\right)={\frac {\partial }{\partial x}}(1)=0\,$

${\frac {\partial ^{2}\theta }{\partial \theta \partial x}}={\frac {\partial ^{2}}{\partial x\partial \theta }}(\theta )={\frac {\partial }{\partial x}}\left({\frac {\partial }{\partial \theta }}(\theta )\right)={\frac {\partial }{\partial x}}(1)=0\,$

${\Big \Downarrow }$

${\frac {\partial ^{2}}{\partial x^{2}}}=\left({\frac {\partial ^{2}}{\partial r^{2}}}{\frac {\partial r}{\partial x}}+{\frac {\partial ^{2}}{\partial r\partial \theta }}{\frac {\partial \theta }{\partial x}}+{\frac {\partial }{\partial \theta }}{\frac {\partial ^{2}\theta }{\partial r\partial x}}\right){\frac {\partial r}{\partial x}}+\left({\frac {\partial ^{2}}{\partial \theta \partial r}}{\frac {\partial r}{\partial x}}+{\frac {\partial }{\partial r}}{\frac {\partial ^{2}r}{\partial \theta \partial x}}+{\frac {\partial ^{2}}{\partial \theta ^{2}}}{\frac {\partial \theta }{\partial x}}\right){\frac {\partial \theta }{\partial x}}\,$

This may be made slightly easier to work with if a few changes are made to the way some of the derivatives are written. Also, the $y$ variable follows analogously:

${\frac {\partial ^{2}}{\partial x^{2}}}=\left({\frac {\partial ^{2}}{\partial r^{2}}}{\frac {\partial r}{\partial x}}+{\frac {\partial ^{2}}{\partial r\partial \theta }}{\frac {\partial \theta }{\partial x}}+{\frac {\partial }{\partial \theta }}{\frac {\partial }{\partial r}}\left({\frac {\partial \theta }{\partial x}}\right)\right){\frac {\partial r}{\partial x}}+\left({\frac {\partial ^{2}}{\partial \theta \partial r}}{\frac {\partial r}{\partial x}}+{\frac {\partial }{\partial r}}{\frac {\partial }{\partial \theta }}\left({\frac {\partial r}{\partial x}}\right)+{\frac {\partial ^{2}}{\partial \theta ^{2}}}{\frac {\partial \theta }{\partial x}}\right){\frac {\partial \theta }{\partial x}}\,$

${\frac {\partial ^{2}}{\partial y^{2}}}=\left({\frac {\partial ^{2}}{\partial r^{2}}}{\frac {\partial r}{\partial y}}+{\frac {\partial ^{2}}{\partial r\partial \theta }}{\frac {\partial \theta }{\partial y}}+{\frac {\partial }{\partial \theta }}{\frac {\partial }{\partial r}}\left({\frac {\partial \theta }{\partial y}}\right)\right){\frac {\partial r}{\partial y}}+\left({\frac {\partial ^{2}}{\partial \theta \partial r}}{\frac {\partial r}{\partial y}}+{\frac {\partial }{\partial r}}{\frac {\partial }{\partial \theta }}\left({\frac {\partial r}{\partial y}}\right)+{\frac {\partial ^{2}}{\partial \theta ^{2}}}{\frac {\partial \theta }{\partial y}}\right){\frac {\partial \theta }{\partial y}}\,$

Now we need to obtain expressions for some of the derivatives appearing above. The most direct path would use the concept of differentials. If:

$x=r\cos(\theta )\quad ,\quad y=r\sin(\theta )\,$

Then:

$dx=dr\cos(\theta )-r\sin(\theta )d\theta \quad ,\quad dy=dr\sin(\theta )+r\cos(\theta )d\theta \,$

Solving by substitution for $dr$ and $d\theta$ gives:

$dr=\cos(\theta )dx+\sin(\theta )dy\quad ,\quad d\theta =-{\frac {\sin(\theta )}{r}}dx+{\frac {\cos(\theta )}{r}}dy\,$

If $z=z(x,y)$, then the total differential is given as:

$dz={\frac {\partial z}{\partial x}}dx+{\frac {\partial z}{\partial y}}dy\,$

Note that the two previous equations are of this form (recall that $r=r(x,y)$ and $\theta =\theta (x,y)$, just like $z$ above), which means that:

$dr={\frac {\partial r}{\partial x}}dx+{\frac {\partial r}{\partial y}}dy\quad ,\quad d\theta ={\frac {\partial \theta }{\partial x}}dx+{\frac {\partial \theta }{\partial y}}dy\,$

Equating coefficients quickly yields a bunch of derivatives:

${\frac {\partial r}{\partial x}}=\cos(\theta )\quad ,\quad {\frac {\partial r}{\partial y}}=\sin(\theta )\quad ,\quad {\frac {\partial \theta }{\partial x}}=-{\frac {\sin(\theta )}{r}}\quad ,\quad {\frac {\partial \theta }{\partial y}}={\frac {\cos(\theta )}{r}}\,$

There's an easier but more abstract way to obtain the derivatives above that may be overkill but is worth mentioning anyway. The Jacobian of the functions $x(r,\theta )$ and $y(r,\theta )$ is:

${\frac {\partial (x,y)}{\partial (r,\theta )}}={\begin{bmatrix}{\cfrac {\partial x}{\partial r}}&{\cfrac {\partial x}{\partial \theta }}\\{\cfrac {\partial y}{\partial r}}&{\cfrac {\partial y}{\partial \theta }}\\\end{bmatrix}}={\begin{bmatrix}\cos(\theta )&-r\sin(\theta )\\\sin(\theta )&r\cos(\theta )\\\end{bmatrix}}$

Note that the Jacobian is a compact representation of the coefficients of the total derivative; using $\mathbf {z} =\mathbf {z} (\mathbf {x} )$ as an example (bold indicating vectors):

$d\mathbf {z} ={\frac {\partial (z_{1},z_{2},...\ ,z_{n})}{\partial (x_{1},x_{2},...\ ,x_{n})}}\cdot d\mathbf {x}$

So, it follows then that the derivatives that we're interested in may be obtained by inverting the Jacobian matrix:

${\frac {\partial (r,\theta )}{\partial (x,y)}}=\left({\frac {\partial (x,y)}{\partial (r,\theta )}}\right)^{-1}$

${\begin{bmatrix}{\cfrac {\partial r}{\partial x}}&{\cfrac {\partial r}{\partial y}}\\{\cfrac {\partial \theta }{\partial x}}&{\cfrac {\partial \theta }{\partial y}}\\\end{bmatrix}}={\begin{bmatrix}{\cfrac {\partial x}{\partial r}}&{\cfrac {\partial x}{\partial \theta }}\\{\cfrac {\partial y}{\partial r}}&{\cfrac {\partial y}{\partial \theta }}\\\end{bmatrix}}^{-1}$

${\begin{bmatrix}{\cfrac {\partial r}{\partial x}}&{\cfrac {\partial r}{\partial y}}\\{\cfrac {\partial \theta }{\partial x}}&{\cfrac {\partial \theta }{\partial y}}\\\end{bmatrix}}={\begin{bmatrix}\cos(\theta )&-r\sin(\theta )\\\sin(\theta )&r\cos(\theta )\\\end{bmatrix}}^{-1}$

${\begin{bmatrix}{\cfrac {\partial r}{\partial x}}&{\cfrac {\partial r}{\partial y}}\\{\cfrac {\partial \theta }{\partial x}}&{\cfrac {\partial \theta }{\partial y}}\\\end{bmatrix}}={\begin{bmatrix}\cos(\theta )&\sin(\theta )\\-{\cfrac {\sin(\theta )}{r}}&{\cfrac {\cos(\theta )}{r}}\\\end{bmatrix}}$

Though somewhat obscure, this is very convenient and it's just one of the many utilities of the Jacobian matrix. An interesting bit of insight is gained: coordinate changes are senseless unless the Jacobian is invertible everywhere except at isolated points, stated another way the determinant of the Jacobian matrix must be nonzero, otherwise the coordinate change is not one-to-one (note that the determinant will be zero at $r=0$ in this example. An isolated point such as this is not problematic.).

Either path you take, there should now be enough information to evaluate the Cartesian second derivatives. Working on $x$:

${\frac {\partial ^{2}}{\partial x^{2}}}=\left({\frac {\partial ^{2}}{\partial r^{2}}}\cos(\theta )-{\frac {\partial ^{2}}{\partial r\partial \theta }}{\frac {\sin(\theta )}{r}}-{\frac {\partial }{\partial \theta }}{\frac {\partial }{\partial r}}\left({\frac {\sin(\theta )}{r}}\right)\right)\cos(\theta )-\left({\frac {\partial ^{2}}{\partial \theta \partial r}}\cos(\theta )+{\frac {\partial }{\partial r}}{\frac {\partial }{\partial \theta }}\left(\cos(\theta )\right)-{\frac {\partial ^{2}}{\partial \theta ^{2}}}{\frac {\sin(\theta )}{r}}\right){\frac {\sin(\theta )}{r}}\,$

${\frac {\partial ^{2}}{\partial x^{2}}}=\cos(\theta )^{2}{\frac {\partial ^{2}}{\partial r^{2}}}-{\frac {\sin(\theta )\cos(\theta )}{r}}{\frac {\partial ^{2}}{\partial r\partial \theta }}+{\frac {\sin(\theta )\cos(\theta )}{r^{2}}}{\frac {\partial }{\partial \theta }}-{\frac {\sin(\theta )\cos(\theta )}{r}}{\frac {\partial ^{2}}{\partial \theta \partial r}}+{\frac {\sin(\theta )^{2}}{r}}{\frac {\partial }{\partial r}}+{\frac {\sin(\theta )^{2}}{r^{2}}}{\frac {\partial ^{2}}{\partial \theta ^{2}}}\,$

${\frac {\partial ^{2}}{\partial x^{2}}}=\cos(\theta )^{2}{\frac {\partial ^{2}}{\partial r^{2}}}-{\frac {2\sin(\theta )\cos(\theta )}{r}}{\frac {\partial ^{2}}{\partial r\partial \theta }}+{\frac {\sin(\theta )\cos(\theta )}{r^{2}}}{\frac {\partial }{\partial \theta }}+{\frac {\sin(\theta )^{2}}{r}}{\frac {\partial }{\partial r}}+{\frac {\sin(\theta )^{2}}{r^{2}}}{\frac {\partial ^{2}}{\partial \theta ^{2}}}\,$

Proceeding similarly for $y$:

${\frac {\partial ^{2}}{\partial y^{2}}}=\left({\frac {\partial ^{2}}{\partial r^{2}}}\sin(\theta )+{\frac {\partial ^{2}}{\partial r\partial \theta }}{\frac {\cos(\theta )}{r}}+{\frac {\partial }{\partial \theta }}{\frac {\partial }{\partial r}}\left({\frac {\cos(\theta )}{r}}\right)\right)\sin(\theta )+\left({\frac {\partial ^{2}}{\partial \theta \partial r}}\sin(\theta )+{\frac {\partial }{\partial r}}{\frac {\partial }{\partial \theta }}\left(\sin(\theta )\right)+{\frac {\partial ^{2}}{\partial \theta ^{2}}}{\frac {\cos(\theta )}{r}}\right){\frac {\cos(\theta )}{r}}\,$

${\frac {\partial ^{2}}{\partial y^{2}}}=\sin(\theta )^{2}{\frac {\partial ^{2}}{\partial r^{2}}}+{\frac {\sin(\theta )\cos(\theta )}{r}}{\frac {\partial ^{2}}{\partial r\partial \theta }}-{\frac {\sin(\theta )\cos(\theta )}{r^{2}}}{\frac {\partial }{\partial \theta }}+{\frac {\sin(\theta )\cos(\theta )}{r}}{\frac {\partial ^{2}}{\partial \theta \partial r}}+{\frac {\cos(\theta )^{2}}{r}}{\frac {\partial }{\partial r}}+{\frac {\cos(\theta )^{2}}{r^{2}}}{\frac {\partial ^{2}}{\partial \theta ^{2}}}\,$

${\frac {\partial ^{2}}{\partial y^{2}}}=\sin(\theta )^{2}{\frac {\partial ^{2}}{\partial r^{2}}}+{\frac {2\sin(\theta )\cos(\theta )}{r}}{\frac {\partial ^{2}}{\partial r\partial \theta }}-{\frac {\sin(\theta )\cos(\theta )}{r^{2}}}{\frac {\partial }{\partial \theta }}+{\frac {\cos(\theta )^{2}}{r}}{\frac {\partial }{\partial r}}+{\frac {\cos(\theta )^{2}}{r^{2}}}{\frac {\partial ^{2}}{\partial \theta ^{2}}}\,$

Now, add these tirelessly hand crafted differential operators and watch the result collapse into just 3 nontrigonometric terms:

$\nabla ^{2}={\frac {\partial ^{2}}{\partial x^{2}}}+{\frac {\partial ^{2}}{\partial y^{2}}}\,$

${\mathrm{\nabla <!--\; \nabla \; -->}}^2=\cos <!--\; \; -->(\mathrm{\theta <!--\; \theta \; -->}{)}^2\frac{{\mathrm{\partial <!--\; \partial \; -->}}^2}{\mathrm{\partial <!--\; \partial \; -->}{r}^2}-<!--\; -\; -->\frac{2\sin <!--\; \; -->(\mathrm{\theta <!--\; \theta \; -->})\cos <!--\; \; -->(\mathrm{\theta <!--\; \theta \; -->})}{r}\frac{{\mathrm{\partial <!--\; \partial \; -->}}^2}{\mathrm{\partial <!--\; \partial \; -->}r\mathrm{\partial <!--\; \partial \; -->}\mathrm{\theta <!--\; \theta \; -->}}+\frac{\sin <!--\; \; -->(\mathrm{\theta <!--\; \theta \; -->})\cos <!--\; \; -->(\mathrm{\theta <!--\; \theta \; -->})}{r^2}\frac{\mathrm{\partial <!--\; \partial \; -->}}{\mathrm{\partial <!--\; \partial \; -->}\mathrm{\theta <!--\; \theta \; -->}}+\frac{\sin <!--\; \; -->(\mathrm{\theta <!--\; \theta \; -->}{)}^2}{r}\frac{\mathrm{\partial <!--\; \partial \; -->}}{\mathrm{\partial <!--\; \partial \; -->}r}+\frac{\sin <!--\; \; -->(\mathrm{\theta <!--\; -->}}{}$