Heat equation with radial symmetry

This notebook contains code to examine Bessel Functions, their orthogonality relation, and their usefulness for solving heat conduction in a plate with radial symmetry.

Bessel Functions

I'm tired of typing this all the time, so let's define a "fat" plotstyle.

fat = AbsoluteThickness[3]

Compare Bessel function (Subscript[J, 0])[r] with a shifted Cos[r]:

Plot[{BesselJ[0,r], Cos[r-Pi/4]},{r,0,25}, PlotStyle -> fat]

Note that if we include the asypmtotic envelope, (Sqrt[2/(πr)]), the functions are visually indistinguishable over the region of interest.

Plot[{BesselJ[0,r], Cos[r-Pi/4]*Sqrt[2/(Pi*r)]},{r,0,25}, PlotStyle -> fat]

Since we will be using the roots of this Bessel function, let's define a shorthand notation for the built-in function BesselJZero[0,n].

k[n_] :=N[BesselJZero[0,n] ]

Table[ k[i], {i,1,3}]

Plot the first 3 eigenfunctions:

Plot[{BesselJ[0, k[1]*r],
BesselJ[0, k[2]*r],
BesselJ[0, k[3]*r]},{r,0,1},
PlotStyle -> fat]

Note that each eigenfunction has one more zero than its predecessor.

Test Orthogonality:

NIntegrate[BesselJ[0,k[1]*r]*BesselJ[0,k[2]*r], {r,0,1}]

That's not zero! What went wrong?
Need weighting factor r:

Integrate[BesselJ[0,k[1]*r]*BesselJ[0,k[2]*r]*r, {r,0,1}]

This is an approximation to zero which you can Chop as desired to make orthogonality more apparent.

Integrate[BesselJ[0,k[1]*r]*BesselJ[0,k[2]*r]*r, {r,0,1}] // Chop

Now let's look at satisfying the initial radial temperature distribution.

Define the initial temperature distribution:

f[r_] := Piecewise[{{1,r<1/2},{ 0,r>1/2}}]

Test that the function produces numerical values before trying to plot:

f[0.2]

Looks good, so now plot initial radial temperature:

Plot[f[r],{r,0,1},PlotStyle -> fat]

Now use the "magic" of projection onto an orthogonal basis to produce the coefficients in the series: (Subscript[c, n])=((f[r]·(Subscript[J, 0])[(Subscript[k, n]) r])/((Subscript[J, 0])[(Subscript[k, n]) r]·(Subscript[J, 0])[(Subscript[k, n]) r]))=((∫f[r] (Subscript[J, 0])[(Subscript[k, n]) r] rdr)/(∫((Subscript[J, 0])^2)[(Subscript[k, n]) r] rdr))

c[n_] := 
Integrate[f[r]*BesselJ[0,k[n]*r]*r, {r,0,1}] /Integrate[BesselJ[0,k[n]*r]*BesselJ[0,k[n]*r]*r, {r,0,1}]

Evaluate the first few coefficients to make sure that things are working OK:

coeffs3 = Table[c[n],{n,1,3}]

Now compute (and save) first 9 coefficients:

coeffs = Table[c[n],{n,1,9}]

Create a corresponding list of orthogonal eigenfunctions to use as basis functions:

basisFns = Table[BesselJ[0,k[n]*r],{n,1,9}]

Plot[basisFns,{r,0,1}]

To check that these can produce a linear combination to approximate the initial conditions, create the linear combination (via inner produce of the lists of coefficients and basis functions), then plot and compare with ICs:

coeffs.basisFns

Plot[{%,f[r]}, {r,0,1}, PlotStyle -> fat]

While not perfect, we already have a reasonable approximation of the ICs. Experiment to see what happens as you keep more terms in the series...

To obtain the time history, we need to include the exponential decay factor that arises from the t-equation. Create a list of "decayFns" that include the exponential time decay along with its associated spatial eigenfunction:

decayFns = Table[BesselJ[0,k[n]*r]*Exp[-k[n]^2*t],{n,1,9}]

Create the series form of the approximate solution by forming a linear combination with the coefficients obtained by projecting onto the orthogonal basis:

soln = coeffs.decayFns

Plot the results for inspection:

Plot3D[soln,{r,0,1},{t,0,0.2},
PlotRange -> {0,1}, AxesLabel -> {x,t,u}]

Think about ways to show the polar geometry of the problem and to better illustrate the evolution. Here is one of many possible approaches:

Animate[
RevolutionPlot3D[
soln/. t->u, 
{r,0,1},
PlotRange -> {-0.1,1.1}],
{u,0,0.025,0.0025}]