Properties of Bessel Functions

Clear["`*"]

Look at the first three Bessel functions as functions of ρ

Plot[{BesselJ[0, x], BesselJ[1, x], BesselJ[2, x]}, {x, 0, 20}, 
      PlotStyle ->  {Hue[0], Hue[0.25], Hue[0.6]}]

Starting with approximate valuesgrom plot, determine accurate values of roots of (Subscript[J, 0])(x),  (Subscript[J, 1])(x),  and (Subscript[J, 2])(x).

x0 = {2.5, 5.5, 8.5, 11.7, 15, 18}

ll  =  Length[x0]

FindRoot[BesselJ[0, x], {x, x0}]

x0 =  x /.  %

x1 = {3.9, 7.0, 10.1, 13.2, 16.3, 19.5}

FindRoot[BesselJ[1, x], {x, x1}]

x1 = x /.  %

x2 =  {5, 8.5, 11.5, 15, 18, 21}

FindRoot[BesselJ[2, x], {x, x2}]

x2  =  x /.  %

Verify that these are indeed roots.

BesselJ[0, x0]

BesselJ[1, x1]

BesselJ[2, x2]

Show that the functions (Subscript[J, n])((Subscript[x, nk]) ρ) and (Subscript[J, n])((Subscript[x, nl]) ρ) are orthogonal with weight ρ on interval [0,1]. Show that the norm is (((1/2)[(Subscript[J, n + 1])((Subscript[x, nk]))])^2).

Table[
        Integrate[
          BesselJ[0, x0[[i]] ρ] BesselJ[0, 
              x0[[j]] ρ] ρ, {ρ, 0, 1}], {i, 1, ll}, {j, 
          1, ll}] // TableForm

Table[BesselJ[1, x0[[i]]]^2/2, {i, 1, ll}]

Table[
        Integrate[
          BesselJ[1, x1[[i]] ρ] BesselJ[1, 
              x1[[j]] ρ] ρ, {ρ, 0, 1}], {i, 1, ll}, {j, 
          1, ll}] // TableForm

Table[BesselJ[2, x1[[i]]]^2/2, {i, 1, ll}]

Table[
        Integrate[
          BesselJ[2, x2[[i]] ρ] BesselJ[2, 
              x2[[j]] ρ] ρ, {ρ, 0, 1}], {i, 1, ll}, {j, 
          1, ll}] // TableForm

Table[BesselJ[3, x2[[i]]]^2/2, {i, 1, ll}]

Plot the functions (Subscript[J, 0])((Subscript[x, 0  k]) ρ), k=1, 2, ... on the interval [0,1]

Plot[{BesselJ[0, x0[[1]] ρ], BesselJ[0, x0[[2]] ρ], 
        BesselJ[0, x0[[3]] ρ], BesselJ[0, x0[[4]] ρ], 
        BesselJ[0, x0[[5]] ρ], 
        BesselJ[0, x0[[6]] ρ]}, {ρ, 0, 1}, 
      PlotStyle ->  {Hue[0], Hue[0.2], Hue[0.4], Hue[0.6], Hue[0.8], 
          Hue[1.0]}, PlotRange -> {-0.5, 1.0}]

Plot the functions (Subscript[J, 1])((Subscript[x, 1  k]) ρ), k=1, 2, ... on the interval [0,1]

Plot[{BesselJ[1, x1[[1]] ρ], BesselJ[1, x1[[2]] ρ], 
        BesselJ[1, x1[[3]] ρ], BesselJ[1, x1[[4]] ρ], 
        BesselJ[1, x1[[5]] ρ], 
        BesselJ[1, x1[[6]] ρ]}, {ρ, 0, 1}, 
      PlotStyle ->  {Hue[0], Hue[0.2], Hue[0.4], Hue[0.6], Hue[0.8], 
          Hue[1.0]}, PlotRange -> {-0.5, 1.0}]

Plot the functions (Subscript[J, 2])((Subscript[x, 2  k]) ρ), k=1, 2, ... on the interval [0,1]

Plot[{BesselJ[2, x2[[1]] ρ], BesselJ[2, x2[[2]] ρ], 
        BesselJ[2, x2[[3]] ρ], BesselJ[2, x2[[4]] ρ], 
        BesselJ[2, x2[[5]] ρ], 
        BesselJ[2, x2[[6]] ρ]}, {ρ, 0, 1}, 
      PlotStyle ->  {Hue[0], Hue[0.2], Hue[0.4], Hue[0.6], Hue[0.8], 
          Hue[1.0]}, PlotRange -> {-0.5, 1.0}]

Cylindrically Symmetric BVP

Clear["`*"]

Consider the case:  Φ(ρ,ϕ,1) = V ( 1- (ρ^2)),   Φ(1,ϕ,z) = 0, and  Φ(ρ,ϕ,0) = 0.  The solution is  Φ(ρ,ϕ,z)  = (∑_(k = 1)^∞) (Subscript[a, 0  k]) (Subscript[J, 0])((Subscript[x, 0  k])ρ) Sinh[(Subscript[x, 0  k])z]/Sinh[(Subscript[x, 0  k])] where (Subscript[a, 0  k])= (({2/[(Subscript[J, 1])((Subscript[x, 0  k]))])^2)} \(TraditionalForm\`∫\_0\%1\( V(
            1 - ρ\^2)\)\ ρ\ \(\(J\_0\)(\(x\_\(0  
                    k\)\) ρ)\)\ \(\(d\)\(ρ\)\(\ \)\)\). Note the following result:

g[a_]  =  
      Integrate[(ρ - ρ^3) BesselJ[0, a ρ], {ρ, 0, 
          1}]

It follows that the expansion coefficient is (Subscript[a, 0  k])= (({4 (Subscript[J, 2])((Subscript[x, 0  k]))/[(Subscript[x, 0  k])(Subscript[J, 1])((Subscript[x, 0  k]))])^2)}.  

Let us determine values of (Subscript[x, 0  k]) over an extended range.

l0  =  30

x0 := {2.5, 5.5, 8.5, 11.7, 15, 18, 
            Table[n π - π/4, {n, 7, l0}]} // Flatten // 
        N

FindRoot[BesselJ[0, x], {x, x0}]

x0 = x /.  %

Define the array a0[k] of expansion coefficients

a0  :=  
      Table[2 g[x0[[k]]]/( BesselJ[1, x0[[k]]])^2, {k, 1, 
          l0}]

Define the potential function

Φ[ρ_, z_]  :=  
      Sum[ a0[ [k]] BesselJ[0, 
            x0[[k]] ρ] Sinh[x0[[k]] z] /
            Sinh[x0[[k]]], {k, 1, l0}]

Check that the potential looks OK on the top surface

<< Graphics`SurfaceOfRevolution`

SurfaceOfRevolution[Φ[ρ, 
        1], {ρ, 0, 1}]

Look at the charge distribution on the top surface

σu[ρ_]  :=  
      D[Φ[ρ, z], z] /.  z ->  1

SurfaceOfRevolution[σu[ρ], {ρ, 0, 1}]

Total charge on upper face (in units (Subscript[ϵ, 0])/V)

Qu = 
      2 π  NIntegrate[ρ σu[ρ], {ρ, 0, 
            1}]

Look at the charge distribution on the lower surface

σl[ρ_]  :=  - D[Φ[ρ, z], z] /.  
        z ->  0

SurfaceOfRevolution[σl[ρ], {ρ, 0, 1}]

Total charge on lower face

Ql  =  
      2 π NIntegrate[ρ σl[ρ], {ρ, 0, 
            1}]

Qu + Ql

Look at the charge distribution on the cylindrical surface

σs[z_]  :=  
      D[Φ[ρ, z], ρ] /.  ρ ->  1

Plot[σs[z], {z, 0, 1}, 
      PlotStyle ->  {Thickness[0.01], RGBColor[0, 1, 0]},  
      AxesLabel -> {z, "σ(z)"}, 
      TextStyle ->  {FontSize ->  14}]

Total charge on cylindrical surface

Qs  =   2  π NIntegrate[σs[z], {z, 0, 1}]

Total charge on all surfaces of cylinder

Qu + Ql + Qs