The Laplacian and heat

For Advanced Calculus

In class we derived the 2D heat equation:

(Subscript[u, t])=k(∇^2)u,

where (∇^2) is the Laplacian. In Cartesian coordinates:

(∇^2)u = (((∂^2)u)/(∂(x^2)))+(((∂^2)u)/(∂(y^2))).

In this notebook, we’ll explore the heat equation together with Laplace’s equation:

(∇^2)u=0.

The 1D heat equation

Well posed problems in one dimension

In one dimension, the Laplacian is just the second derivative. Thus, the heat equation becomes

(Subscript[u, t])=k (((d^2)u)/(d(x^2))).

Generally, PDEs are not of interest unless unaccompanied by other conditions. In the simplest case, we study the diffusion of heat through a bar of length L lying along the interval [0,L]. The function u(x,t) represents the temperature of the bar at position x and time t. The condition on u that we need, in addition to the heat equation itself, are initial and boundary conditions. The initial condition is specified by a function of position. We write

u(x,0)=f(x).

The boundary conditions are specified as functions of t. We write

u(0,t)=(Subscript[g, 0])(t)  and  u(L,t)=(Subscript[g, L])(t).

Quite often, these might just be constants.

An example

An example of a fully written out heat problem might look like so:

({{(Subscript[u, t])=(Subscript[u, x x]), u(x,0)=27(x^2)(1-x)/4}, {u(0,t)=0, u(1,t)=0}})

Note that k=1, L=1, f(x)=27(x^2)(1-x)/4, (Subscript[g, 0])(t)=0, and (Subscript[g, 1])(t)=0. Thus, the initial temperature distribution in the bar looks like so:

Plot[27(x^2)(1-x)/4,{x,0,1}, 
PlotRange -> {-0.05,1.05},PlotStyle->Thick,
ColorFunction->"TemperatureMap",
ColorFunctionScaling->False,
Epilog -> {PointSize[Large],
{ColorData["TemperatureMap"][0],Point[{{0,0},{1,0}}]}
}]

The function u(x,t) tells us how that temperature profile evolves with respect to time. Since the ends are held at constant temperature zero, basic physical principles tell us that we expect the heat to flow out the ends of the bar as the system approaches a steady state of constant zero temperature.

Visualization

We can set up and solve the system from the previous section using a numerical PDE solver. Here’s the system again:

({{(Subscript[u, t])=(Subscript[u, x x]), u(x,0)=27(x^2)(1-x)/4}, {u(0,t)=0, u(1,t)=0}})

Here’s a representation of that same system in Mathematica speak:

Clear[u,heatSystem];
heatSystem = {
D[u[x,t],t]==D[u[x,t],x,x],
u[x,0]==27x^2(1-x)/4,
u[0,t]==0, u[1,t]==0};

We can solver this with Mathematica’s awesome numerical PDE solver:

Clear[u];
u[x_,t_] = NDSolveValue[heatSystem,
u[x,t],{x,0,1},{t,0,1}]

The result is an interpolating function. While it is a black box, we can still compute numerical values with it. Here’s the temperature at the midpoint of the bar after 1 second:

u[0.5,1]

Very close to zero. This makes sense, if we think in terms of steady state.

We can animate the evolution of the system using the groovy Manipulate command.

Manipulate[Plot[u[x,t],{x,0,1},
ColorFunction->"TemperatureMap", 
ColorFunctionScaling->False,
PlotRange -> {0,1}],
{t,0,1/2}]

A little variant on this creates an animated GIF that you can upload to our discussion forum:

pics = Table[Plot[u[x,t],{x,0,1},
ColorFunction->"TemperatureMap", 
ColorFunctionScaling->False,
PlotRange -> {0,1}],
{t,0,0.5,0.01}];
Export["anim.gif",pics]

Consistency

In the previous example, we had u(x,0)=27(x^2)(1-x). Thus, u(0,0)=u(1,0)=0, in agreement with the boundary conditions.  We say that the initial condition and the boundary conditions are consistent. Here’s an example of an inconsistent system:

({{(Subscript[u, t])=(Subscript[u, x x]), u(x,0)=0}, {u(0,t)=0, u(1,t)=1}})

Note that at the right hand endpoint, the initial condition is u(0,0)=0 but the boundary condition at time zero is u(1,0)=1. Thus the initial and boundary conditions are inconsistent. This situation does arise in practice. For example, one end of a hot bar might be brought into immediate contact with a hot stove. Mathematica’s numerical PDE solver can deal with this, though some special care is required.

Clear[u];
u[x_,t_] = NDSolveValue[{
D[u[x,t],t]==D[u[x,t],{x,2}],u[x,0]==0,
u[0,t]==0,u[1,t]==1},
u[x,t],{x,0,1},{t,0,1},
Method->{"MethodOfLines", 
"DifferentiateBoundaryConditions"->{True, "ScaleFactor"->100}}];
Manipulate[Plot[u[x,t],{x,0,1},
ColorFunction->"TemperatureMap", 
ColorFunctionScaling->False,
PlotRange -> {0,1}],
{t,0,1}]

An ask and answer exercise - invite only!!

Use Mathematica to find a numerical solution to the following system:

({{(Subscript[u, t])=(Subscript[u, x x]), u(x,0)=sin(5x)}, {u(0,t)=0, u(1,t)=1}})

Laplace’s equation and steady state

The example in the previous section illustrates a common occurrence  for the heat equation with constant boundary conditions - namely, the system approaches a steady state where the temperature is constant with respect to time. Thus, in the limit, we have (Subscript[u, t])=0 so that the heat equation

(Subscript[u, t])=k(∇^2)u,

becomes Laplace’s equation:

(∇^2)u=0.

Laplace’s equation is an important PDE in its own right. Among other things, solutions to Laplace’s equation with appropriate boundary conditions represent steady state heat flow.

An ask and answer exercise - Invite only!!

Find the solution to the following system consisting of Laplace’s equation with boundary conditions:

(∇^2)u=0, u(0)=-1,u(2)=3.

Your solution should be quite a simple algebraic expression. Graph the solution and interpret it in the context of steady state heat flow.

The polar Laplacian

In polar coordinates, the Laplacian has the form:

(∇^2)U = (Subscript[U, r r])+(1/r)(Subscript[U, r])+(1/(r^2))(Subscript[U, θθ])

Mathematica “knows” this.

Laplacian[U[r,t],{r,t},"Polar"]//Expand

Changing coordinates

By hand

To explain the polar Laplacian, we can carry out a change of coordinates from Cartesian to polar.  Suppose that u is a function of the Cartesian variables x and y and U is a function of the polar variables r and θ related to u by

U(r,θ)=u(r cos(θ),r sin(θ)).

We can use the multi-variable chain rule to compute the polar derivatives of U in terms of the Cartesian derivatives of u.  Recall that the multi-variable chain rule states that if f=f(x,y) is a function of x and y, and x=x(t) and y=y(t) are functions of t, then

(d/(dt))f(x(t),y(t)) = ((∂f)/(∂x))((dx)/(dt))+((∂f)/(∂y))((dy)/(dt)).

We can apply this to U to get the following values of its derivatives.

(Subscript[U, r])  = (Subscript[u, x])cos(θ)+(Subscript[u, y])sin(θ)
(Subscript[U, r r])=((Subscript[u, x x])cos(θ)+(Subscript[u, x y])sin(θ))cos(θ) + ((Subscript[u, y x])cos(θ)+(Subscript[u, y y])sin(θ))sin(θ)
(Subscript[U, θ])  = (Subscript[u, x])(-r sin(θ))+(Subscript[u, y])r cos(θ)
(Subscript[U, θθ])=((Subscript[u, x x])(-r sin(θ))+(Subscript[u, x y])r cos(θ))(-r sin(θ)) + (Subscript[u, x])(-r cos(θ)) +
((Subscript[u, y x])(-r sin(θ))+(Subscript[u, y y])(r cos(θ)))r cos(θ) + (Subscript[u, y])(-r sin(θ))

We now use these to compute (Subscript[U, r r])+(1/r)(Subscript[U, r])+(1/(r^2))(Subscript[U, θθ]) in terms of u.

(Subscript[U, r r])+(1/r)(Subscript[U, r])+(1/(r^2))(Subscript[U, θθ]) = 
(Subscript[u, x x])(cos^2)(θ)+(Subscript[u, x y])sin(θ)cos(θ)+(Subscript[u, y x])cos(θ)sin(θ)+(Subscript[u, y y])(sin^2)(θ) +
(1/r)((Subscript[u, x])cos(θ)+(Subscript[u, y])sin(θ)) +
(1/(r^2))((Subscript[u, x x])(r^2)(sin^2)(θ)-(Subscript[u, x y])(r^2)cos(θ)sin(θ) - (Subscript[u, x])r cos(θ) -
(Subscript[u, y x])(r^2)sin(θ)cos(θ)+(Subscript[u, y y])(r^2)(cos^2)(θ) - (Subscript[u, y])r sin(θ))=
(Subscript[u, x x])((cos^2)(θ)+(sin^2)(θ)) + (Subscript[u, y y])((sin^2)(θ)+(cos^2)(θ)) = (Subscript[u, x x])+(Subscript[u, y y])

In particular,  (Subscript[U, r r])+(1/r)(Subscript[U, r])+(1/(r^2))(Subscript[U, θθ])=0, if (Subscript[u, x x])+(Subscript[u, y y])=0.

Automated

Here's an automated approach.

Clear[u,U];
U[r_,θ_] := u[r*Cos[θ],r*Sin[θ]];
Simplify[(∂_(r,r))U[r,θ]+(1/r)(∂_r)U[r,θ]+(1/(r^2))(∂_(θ,θ))U[r,θ]] /. 
{r*Cos[θ]->x,r*Sin[θ]->y}

We can even go the other way - which is much harder.

Clear[u,U];
u[x_,y_] := U[Sqrt[x^2+y^2],ArcTan[y/x]];
(FullSimplify[(∂_(x,x))u[x,y]+(∂_(y,y))u[x,y]] /. 
{(x^2)+(y^2)->(r^2),ArcTan[y/x]->θ}) /.{((((r^2)))^(1/2)) -> r,((((r^2)))^(-1/2))  -> 1/r}

A simple example

The function u(x,y)=4(x^3)y-4x (y^3) is harmonic on the plane, as you can easily verify.

u[x_,y_] := 4(x^3)y-4x (y^3);
(∂_(x,x))u[x,y] + (∂_(y,y))u[x,y]

Here is the same function expressed in polar coordinates.

U[r_,θ_] = u[r Cos[θ],r Sin[θ]]

There is no reason to expect that (Subscript[U, r r])+(Subscript[U, θθ])=0.

(∂_(r,r))U[r,θ]+(∂_(θ,θ))U[r,θ]

But we do get 0, if we apply the polar Laplacian.

(∂_(r,r))U[r,θ]+(1/r)(∂_r)U[r,θ]+(1/(r^2))(∂_(θ,θ))U[r,θ]//Expand

Now that we have done this, we will continue to typically denote the state variable of a PDE by a lower-case u, regardless of the coordinate system used.

Radial symmetry

When u has rotational symmetry, its Laplacian simplifies to

(∇^2)u = (Subscript[u, r r])+(1/r)(Subscript[u, r]).

This is an ODE, that can be reduced to first order with the substitution y=(Subscript[u, r]).  It’s not so bad to find it’s solutions.

Clear[u];
DSolve[u''[r]+u'[r]/r==0,u[r],r]

Here’s a specific solution on an annulus, for example.

Clear[u];
DSolve[{u''[r]+u'[r]/r==0,u[1]==0,u[2]==2},u[r],r]

Let’s visualize.

ParametricPlot3D[{r*Cos[θ],r*Sin[θ],2Log[r]/Log[2]},{r,1,2},{θ,0,2π},
Mesh -> False,BoundaryStyle->Thick,PlotPoints -> 100,
ColorFunctionScaling->False,ColorFunction -> "TemperatureMap"]

Or

ParametricPlot[{r*Cos[θ],r*Sin[θ]},{r,1,2},{θ,0,2π},
Axes -> False,BoundaryStyle -> None,PlotPoints -> 100,
Mesh -> False,ColorFunctionScaling->False,ColorFunction -> 
Function[{x,y,r},ColorData["TemperatureMap"][2Log[r]/Log[2]]]]

Or

DensityPlot[Log[(x^2)+(y^2)]/Log[2],{x,-2.2,2.2},{y,-2.2,2.2},
RegionFunction -> Function[{x,y},1≤x^2+y^2≤4],
ColorFunctionScaling->False,
ColorFunction->"TemperatureMap", PlotPoints -> 100]

The 2D heat equation

Cartesian coordinates

The 2D heat equation is similar to the 1D but requires more boundary conditions. The following code specifies that the solution should be zero on the boundary of the unit square. The initial condition is a bit rigged to be consistent with the boundary conditions.

Clear[u];
u[x_,y_,t_] = NDSolveValue[{
D[u[x,y,t],t]==Laplacian[u[x,y,t],{x,y}],
u[x,0,t]==0,u[1,y,t]==0,
u[x,1,t]==0,u[0,y,t]==0,
u[x,y,0] == 16x(1-x)y(1-y)
},u[x,y,t],{x,0,1},{y,0,1}, {t,0,1/2}]

We can visualize the solution as a 3D plot:

Manipulate[Plot3D[u[x,y,t],{x,0,1},{y,0,1},
ColorFunction->"TemperatureMap",
ColorFunctionScaling->False, PerformanceGoal->"Quality",
PlotRange -> {0,1},
MeshFunctions->{#3&}, Mesh -> {Range[0,1,0.1]}],
{t,0,0.2}]

Or as a density plot:

Manipulate[DensityPlot[u[x,y,t],{x,0,1},{y,0,1},
ColorFunction->"TemperatureMap",
ColorFunctionScaling->False, PerformanceGoal->"Quality",
PlotPoints -> 32,
MeshFunctions->{#3&}, Mesh -> {Range[0,1,0.1]}],
{t,0,0.2}]

Here’s a similar example with inconsistent boundary conditions because one side is changed instantaneously.

Clear[u];
u[x_,y_,t_] = NDSolveValue[{
D[u[x,y,t],t]==Laplacian[u[x,y,t],{x,y}],
u[x,0,t]==0,u[1,y,t]==1,
u[x,1,t]==0,u[0,y,t]==0,
u[x,y,0] == 0
},u[x,y,t],{x,0,1},{y,0,1}, {t,0,1/2},
Method->{"MethodOfLines", 
"DifferentiateBoundaryConditions"->{True, "ScaleFactor"->100}}]

Polar coordinates

The pringle

You can check that u(r,θ)=(r^2)sin(2θ) satisfies (∇^2)u=0. Here’s what it looks like:

boundary = ParametricPlot3D[{Cos[θ],Sin[θ],Sin[2θ]},
{r,0,1},{θ,0,2π}, PlotStyle -> {Black,Thickness[0.02]}];
interior = ParametricPlot3D[{r*Cos[θ],r*Sin[θ],(r^2)Sin[2θ]},
{r,0,1},{θ,0,2π}, BoundaryStyle -> Black, 
ColorFunction -> "TemperatureMap",MeshFunctions -> {#3&}];
Show[{interior, boundary}]

The Poisson integral formula

The Poisson integral formula states that

u(r,θ)=(1/(2π))(∫_(-π)^π)f(ϕ)(((R^2)-(r^2))/((R^2)+(r^2)-2R r cos(θ-ϕ)))dϕ

satisfies Laplace’s equation with u(R,θ)=f(θ) on the boundary of the disk of radius R. We can check that the integrand is harmonic by straight computation.

(Subscript[P, R_])[r_,θ_] = (1/(2π))((((R^2)-(r^2)))/((R^2)+(r^2)-2R r Cos[θ]));
Laplacian[(Subscript[P, R])[r,θ],{r,θ},"Polar"]//FullSimplify

We can apply it like so:

Clear[u];
u[1,θ_] := Sin[2θ];
u[r_,θ_] := NIntegrate[(1/(2π))Sin[2ϕ](((1-(r^2)))/(1+(r^2)-2r Cos[θ-ϕ])),
{ϕ,-π,π}]/; r<1;
pringle =ParametricPlot3D[{r*Cos[θ],r*Sin[θ],u[r,θ]},
{r,0,0.99},{θ,-Pi,Pi}, ColorFunction -> "TemperatureMap",
PlotPoints -> {10,32},MaxRecursion -> 5];
boundary = ParametricPlot3D[0.99{Cos[θ],Sin[θ],Sin[2θ]},{θ,-Pi,Pi},
PlotStyle -> Directive[Thick,Black]];
Show[{pringle,boundary}]

Using NDSolve

Using NDSolve over a disk is a little tricky but quite elegant, once you get the hang of it.

Clear[u];
omega = Disk[{0,0},1];
boundaryCondition = DirichletCondition[u[x,y]==2x*y, x^2+y^2==1];
u[x_,y_] = NDSolveValue[
{Laplacian[u[x,y],{x,y}]==0,boundaryCondition},
u[x,y],Element[{x,y},omega]];

Plot3D[u[x,y],Element[{x,y},omega], Mesh -> All,
ColorFunction -> "TemperatureMap", BoxRatios -> Automatic]

It turns out that you can solve PDEs over quite general regions.

omega = RegionDifference[
Polygon[{{-2,-2},{2,-2},{2,2},{-2,2}}],
Disk[{0,0},1]
];
RegionPlot[omega, PlotRange -> {{-3,3},{-3,3}}]

Let’s set the temperature on the boundary of the interior disk to be 1 and the temperature on the exterior square to be zero.

boundaryConditions = {
DirichletCondition[u[x,y]==1, x^2+y^2==1],
DirichletCondition[u[x,y]==0,x==-2||x==2||y==-2||y==2]
};

Then, we’ll solve the Laplace’s equation over that region with those boundary conditions

Clear[u];
u[x_,y_] = NDSolveValue[
{Laplacian[u[x,y],{x,y}]==0,boundaryConditions},
u[x,y],Element[{x,y},omega]]

And visualize:

Plot3D[u[x,y],Element[{x,y},omega], Mesh -> All,
PlotRange -> {{-2,2},{-2,2},{0,1}},
ColorFunction -> "TemperatureMap", BoxRatios -> Automatic,
BoundaryStyle->{Thick,Black}]

The spherical Laplacian

The algebraic form

Here’s the rather complicated looking spherical Laplacian:

Clear[u];
Laplacian[u[ρ,φ,θ],{ρ,φ,θ},"Spherical"]

I guess we can prove this like so:

Clear[u,U];
U[ρ_,φ_,θ_] := u[ρ*Sin[φ]*Cos[θ],ρ*Sin[φ]*Sin[θ],ρ*Cos[φ]];
Simplify[
Derivative[2,0,0][U][ρ,φ,θ]+((((Derivative[0,2,0][U][ρ,φ,θ])/ρ)+Derivative[1,0,0][U][ρ,φ,θ])/ρ)+((Csc[φ] (((Csc[φ] Derivative[0,0,2][U][ρ,φ,θ])/ρ)+((Cos[φ] Derivative[0,1,0][U][ρ,φ,θ])/ρ)+Sin[φ] Derivative[1,0,0][U][ρ,φ,θ]))/ρ)
] /. 
{ρ*Sin[φ]*Cos[θ]->x,ρ*Sin[φ]*Sin[θ]->y,ρ*Cos[φ]->z}

Like the polar Laplacian, this is much simpler when we have radial symmetry!

Laplacian[u[ρ],{ρ,φ,θ},"Spherical"]

Example

Here’s an example involving Poisson’s equation, (∇^2)u=-k, which describes steady state heat flow with a source.

Clear[u]
DSolve[{u''[r]+2u'[r]/r==-6,u[0]==1},u[r],r]

I guess this says that, if we hold the origin at a constant temperature in a medium with a uniform heat source that diffuses heat uniformly in all directions, then the steady state heat distribution will decrease quadratically as we move away from the origin.