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.