Symbolic solution for steady-state heat equation i.e. Laplace equation inside cylinder

Question

I want to symbolically solve the following boundary value problem in my textbook. It's steady-state heat conduction equation i.e. Laplace equation inside a cylinder

$$ \left\{\begin{array}{l} \Delta u=0\\ \left.u\right|_{\rho=a}=u_{0} \\ \left.u\right|_{z=0}=\left.0 \quad u\right|_{z=h}=0 \end{array}\right. $$

The following series solution can be found manually by method of separation of variables:

$$ \begin{aligned} u&=\sum_{m=1}^{\infty} A_{m} I_{0}\left(\frac{m \pi}{h} \rho\right) \sin \frac{m \pi}{h} z \\ &=\sum_{n=0}^{\infty} \frac{4 u_{0}}{I_{0}\left[\frac{(2 n+1) \pi}{h} a\right](2 n+1) \pi} I_{0}\left[\frac{(2 n+1) \pi}{h} \rho\right] \sin \left[\frac{(2 n+1) \pi}{h} z\right] \end{aligned} $$

But DSolve can't solve the problem:

eq = D[u[rho, z], {rho, 2}] + D[u[rho, z], {rho, 1}]/rho + 
   D[u[rho, z], {z, 2}] == 0
ic1 = {u[rho, 0] == 0, u[rho, 1] == 0};
ic2 = {u[1, z] == 1};

DSolve[{eq, ic1, ic2}, u, {rho, z}]

Why?

Answer 1

This post contains several code blocks, you can copy them easily with the help of functions here.

---

It's not too surprising to see DSolve failing on the problem, because though DSolve is improved these years, it's still fragile. So let me show a solution based on finite Fourier sine transform, which is much cleaner and easier to understand than separation of variables.

I'll use finiteFourierSinTransform for the task.

It's not necessary, but we define 2 functions to make the output pretty:

Format@finiteFourierSinTransform[f_, __] := Subscript[ℱ, s][f]    
Format@u0 := Subscript[u, 0]

Then we write down the equation and transform. Here I'd like to point out that, except for the 3 b.c.s in the question, there actually exists an implicit constraint i.e. the solution is bounded for $\rho<a$. In your specific case, the solution is axisymetric, so the constraint is equivalent to

$$\left.\frac{\partial u}{\partial \rho}\right|_{\rho=0}=0$$

I'll include this new b.c. in the code to facilitate subsequent discussion:

With[{u = u[ρ, z]}, eq = Laplacian[u, {ρ, ϕ, z}, "Cylindrical"] == 0;
  bcz = u == 0 /. {{z -> 0}, {z -> h}};
  bcrexplicit = u == u0 /. ρ -> a;
  bcrimplicit = D[u, ρ] == 0 /. ρ -> 0];

finiteFourierSinTransform[{eq, bcrexplicit, bcrimplicit}, {z, 0, h}, m]

% /. Rule @@@ bcz

tset = % /. HoldPattern@finiteFourierSinTransform[f_, __] :> f /. u -> (U@# &)

Here I use U to represent finite Fourier sine transform of u to facilitate subsequent coding.

tset is 2nd order ODE with 2 b.c.s, and DSolve can handle it:

tsol = DSolveValue[tset, U@ρ, ρ] // FullSimplify

Remark

If we haven't included bcrimplicit in the code, it's still possible to obtain tsol, but the discussion is a bit involved. Without bcrimplicit, we'll need to solve 2nd order ODE with 1 b.c. i.e.

tsolmid = DSolveValue[tset // Most, U@ρ, ρ] // FullSimplify

There remains a constant C[1], which is expected, because we haven't use the implicit constraint that the solution is bounded for $\rho<a$. Sadly there doesn't seem to be a straightforward way to impose this condition in DSolve at this stage, but by picking up, or simply testing properties of BesselY a bit, we know that

BesselY[0, -((I m π ρ)/h)] /. ρ -> 0
(* -∞ *)

Abs@BesselY[0, -((I a m π)/h)] < Infinity // Simplify
(* True *)

So, to keep the solution bounded for $\rho<a$, the coefficient of BesselY[0, -((I m π ρ)/h)] must be 0 i.e.

const = Solve[-(((-1 + (-1)^m) h u0)/(m π)) - BesselI[0, (a m π)/h] C[1] == 
    0, C[1], Assumptions -> h > 0][[1]]

Substitute it back to tsolmid, we obtain the tsol as shown above:

tsol = tsolmid /. const // Simplify

The final step is to transform back:

sol = inverseFiniteFourierSinTransform[tsol, m, {z, 0, h}]

It's not hard to notice this is already the solution in your question. (Notice I use C to represent Infinity when designing inverseFiniteFourierSinTransform. ) You can make it more pretty, of course:

solpretty = 
 sol /. coef_ HoldForm[Sum[expr_, {i_, C}]] :> 
    Inactive[Sum][expr coef, {i, 1, Infinity}] // TraditionalForm

solfinal = solpretty /. m -> 2 n + 1 /. {2 n + 1, 1, ∞} -> {n, 0, ∞} // 
 Simplify[#, n ∈ NonNegativeIntegers] &

(*
solfinal = 
 Inactive[Sum][(
  4 u0 BesselI[0, ((1 + 2 n) π ρ)/h] Sin[((1 + 2 n) π z)/
    h])/((1 + 2 n) π BesselI[0, (a (1 + 2 n) π)/h]), {n, 0, ∞}]
 *)

You can use the solution for visualization, of course:

hvalue = 1; u0value = 1; avalue = 1; nvalue = 10;

asol = solfinal /. {TraditionalForm -> Identity, ∞ -> nvalue, 
     a -> avalue, u0 -> u0value, h -> hvalue} // Activate;

SliceDensityPlot3D[
  asol /. ρ -> Sqrt[x^2 + y^2] // Evaluate, 
  "CenterPlanes", {x, y, z} ∈ Cylinder[{{0, 0, 0}, {0, 0, hvalue}}, avalue]]

Compare it with FEM-based numeric solution:

fem[measure_ : Automatic] := {"FiniteElement", 
        "MeshOptions" -> MaxCellMeasure -> measure};

nsol = NDSolveValue[{eq, bcz, bcrexplicit} /. {a -> avalue, u0 -> u0value, 
     h -> hvalue}, u, {ρ, 0, avalue}, {z, 0, hvalue}, Method -> fem[10^-4]];

Manipulate[
 Plot[{#, nsol[ρ, z]} // Evaluate, {z, 0, hvalue}, 
  PlotStyle -> {Automatic, Dashed}, PlotRange -> {0, 1.5}], {ρ, 0, avalue}] &@asol

The error is a bit obvious at $\rho=a$, but as an approximation with the first 10 terms of the series only, this isn't bad.