3
$\begingroup$

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}]

enter image description here

Why?

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

1 Answer 1

7
$\begingroup$

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]

enter image description here

% /. Rule @@@ bcz

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here

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]]

enter image description here

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}]

enter image description here

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

enter image description here

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, ∞}]
 *)

enter image description here

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]]

enter image description here

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

enter image description here

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.