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

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?

• 2nd $∑_{m=1}^∞$ should be $∑_{n=0}^∞$ , right? Jun 25 at 12:53
• @xzczd m=2n+1,You can do it by ha nd, but you can't do it by software Jun 25 at 13:17
• I mean, there's a typo in your $\LaTeX$, 2nd $m=1$ should be $n=0$. Please correct it. Jun 25 at 13:22
• @xzczd I have correct it Jun 25 at 13:28

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. 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. 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, 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.

• you say finiteFourierSinTransform,This is the first time I've heard of it. I'm going to learn it Jun 25 at 13:31
• @我心永恒 Since you can read Chinese, you may want to read this: zhihu.com/question/412203550/answer/1388426673 Jun 25 at 13:38
• thanks you,I will read it. Jun 25 at 13:41