Generalizations of the Fundamental Theorem of Calculus (Part II)

Using Divergence theorem to find London attraction potential between point particle and a finite cylinder.

\(The\ \(dipole--\) 
      dipole\ interaction\ force\ goes\ like\ \(1\/r\^6\) . \ \ The\ London\ interaction\ that\ gives\ the\ \ force\ between\ two\ bodies\ due\ to\ the\ superposition\ of\ all\ pointwise\ dipole\ interactions . \ 

To\ find\ such\ forces, \ 
    it\ is\ necessary\ to\ integrate\ over\ the\ volumes\ of\ both\ \(bodies--\) an\ integral\ over\ six\ dimensions . \ \ There\ are\ a\ few\ \(cases--\) such\ as\ two\ spheres, \ 
    a\ sphere\ and\ a\ semi - infinite\ solid, \ 
    etc . \(--\(-where\)\)\ closed\ form\ solutions\ can\ be\ obtained . \ \ \ However, \ in\ general, \ 
    it\ is\ necessary\ to\ perform\ a\ numerical\ integration . \ \ If\ the\ integral\ over\ both\ volumes\ could\ be\ converted\ to\ integrals\ over\ both\ surfaces, 
    then\ the\ numerical\ integration\ is\ considerably\ more\ efficient . \ \ Such\ conversions, \ 
    can\ be\ obtained\ though\ application\ of\ the\ divergence\ theorem . 

\ 
        Here\ is\ an\ example\ for\ the\ interaction\ between\ a\ single\ \((point)\)\ and\ a\ finite\ volume\ \((a\ cylinder)\) . \ \ To\ proceed, \ \ a\ vector\ \ integrand\ \ must\ be\ found\ for\ which\ the\ divergence\ is\ the\ desired\ \ pointwise\ \(\(interaction\)\(.\)\)\)

Find (Overscript[F,  ->])  such that div (Overscript[F,  ->])  is \(TraditionalForm\`\(-1\)\/\((\(r\& -> \) - \(x\& -> \))\)\^6\)   where (Overscript[r,  ->])  = (ξ,η,ζ) is a position in the cylinder and (Overscript[x,  ->]) =(x,y,z) is a general position in space (A fairly general method to do this and similar problems can be found in Argento C; Jagota A; Carter WC ``Surface formulation for molecular interactions of macroscopic bodies'' J.Mech. Physics Solids 1997,  pp 1161-1183 .

The following is a ``guess'' at the vector potential; it will be verified as the correct one by checking its divergence.

FVecLondon  =  (1/(3 ((((((ξ - x))^2)  + (((η - y))^2)  +  (((ζ - z))^2)))^3))) {ξ - x,  η - y,  ζ - 
            z}

<< Calculus`VectorAnalysis`

The following verifies that the correct vector potential is obtained

FullSimplify[
      Div[FVecLondon, Cartesian[ξ, η, ζ]]]

Therefore,  \(TraditionalForm\`∫\_\(Cylinder

          Volume\)\(1\/\((\(r\& -> \) - \(x\& -> \))\)\^6\) dV\) = \(TraditionalForm\`∫\_\(Cylinder

          Surfaces\)\(FVecLondon\&⇀\) • d\(A\& -> \)\) =  \(TraditionalForm\`∫\_\(Cylinder

          Surface\)\(FVecLondon\&⇀\) • d\(A\& -> \)\)+ \(TraditionalForm\`∫\_\(Cylinder

          Ends\)\(FVecLondon\&⇀\) • d\(A\& -> \)\) is the total interaction between a point an a cylinder.

Cylindrical coordinates: (ξ,η,ζ) = (r Cos[t], r Sin[t], ζ):
One cylindrical surface  r= R = const. t ∈ (0, 2π), ζ ∈ (-(L/2) , (L/2))
Two caps r ∈ (0,R) t ∈ (0, 2π), ζ=±(L/2)

Cylinder Surface normals and differential quantities

The following is a parametric representation of a cylinder surface that is coaxial with the z-axis (the cylinder ends will be included later)

CylSurf  =  {R Cos[t],  R  Sin[t], ζ}

The infinitessimal surface vectors (Overscript[((Subscript[R, u])), ⇀]) and (Overscript[((Subscript[R, v])), ⇀]) for the cylinder surface are obtained by differentiation; they will be used to find the surface patch d(Overscript[A,  ->]) .

CylSurfRt  =  D[CylSurf, t]
CylSurfRz    =  D[CylSurf, ζ]

The surface normal given by (Overscript[((Subscript[R, u])), ⇀]) × (Overscript[((Subscript[R, v])), ⇀]) for the cylinder surface, there for the following (multiplied by dθ dz) is the infinitessimal oriented surface patch d(Overscript[A,  ->]) .

NormalVecCylSurf  =  Cross[CylSurfRt, CylSurfRz]

The integrand to be evaluated over the cylinder surface is the vector potential, dotted into

CylinderIntegrandθζ  =  
      FullSimplify[(FVecLondon /. {ξ -> 
                  R Cos[t],  η ->  R Sin[t]}) . 
          NormalVecCylSurf]

Closed form for integral is un-findable, use a numeric integration--this will be written in two forms
1) A general function of x,y,z of the point particle and the cylinder radius and length.
2) Because of the polar symmetry, the integral should  be a function of only r = (Sqrt[(x^2)  +  (y^2)])

We can use the same function name twice, but with a different number of arguments.  But, first we will try to do at least one of the two integrations symbolically.  If we can do this, it will reduce the original numerical integration from three dimensions down to 1---which is a nice bargain if we can find it.

Try and see if we can do one of the two integrals.  First, let's try theta:

Integrate[CylinderIntegrandθζ, t,  
      Assumptions -> 
        R  >  0  &&   ζ  ∈  Reals  &&  
          x  ∈  Reals && y ∈  Reals &&  
          z  ∈  Reals]

The integral over theta gives no symbolic form...Try and see if we can do one of the two integrals.  So, let's try ζ:

CylinderIntegrandθIndζ = 
      Integrate[CylinderIntegrandθζ, ζ ,  
        Assumptions -> 
          R  >  0  &&   L > 0  &&   x  ∈  Reals && 
            y ∈  Reals &&  z  ∈  Reals &&  
            t ∈  Reals]

Here we will use a trickof using Evaluate[] in a function definition to save time--a small diversion here will demonstrate why this is efficient

Timing[
      f[c_]  :=  Evaluate[Integrate[Exp[Tan[x]], {x, 0, c}]]]

Timing[g[c_]  :=  Integrate[Exp[Tan[x]], {x, 0, c}]]

?f

?g

Timing[f[0.5]]

Timing[g[0.5]]

CylinderIntegrandθ[x_, y_, z_, CylRad_, CylLen_]  :=  
      Evaluate[Simplify[(CylinderIntegrandθIndζ /. {ζ -> 
                    CylLen/2, 
                  R -> CylRad}) - (CylinderIntegrandθIndζ /. {ζ -> -CylLen/2, R -> CylRad}), 
          Assumptions -> 
            CylRad  >  0  &&   CylLen > 0  &&   x  ∈  
                Reals && y ∈  Reals &&  
              z  ∈  Reals &&  t ∈  Reals]]

CylinderContribution[x_,  y_,  z_,  CylRad_ ,  CylLen_]  :=  
      NIntegrate[
        Evaluate[
          CylinderIntegrandθ[x, y, z, CylRad, CylLen], {t, 0,  
            2  π}]]

\(Because\ of\ the\ polar\ symmetry, \ 
    this\ contribution\ should\ only\ be\ a\ function\ of\ the\ distance\ \((\@\(x\^2 + y\^2\))\), \ 
    so\ we\ can\ write\ a\ new\ function\ \((with\ the\ same\ name)\)\ with\ fewer\ arguments\)

CylinderContribution[dist_,  z_,  CylRad_ ,  CylLen_]  := 
      NIntegrate[
        Evaluate[
          CylinderIntegrandθ[dist, 0, z, CylRad, CylLen], {t, 0,  
            2  π}]]

Checking the numerical integration for the first form of the function

CylinderContribution[2, 2, 1/2, 1, 1]

Checking the numerical integration for the second form of the function

CylinderContribution[2 Sqrt[2], 1/2, 1, 1]

Top Cap Surface normals and differential quantities, parallel method to cylinder surface

TopSurf  =  {r Cos[t], r Sin[t], (L/2)}

TopSurfRt  =  D[TopSurf, t]
TopSurfRr    =  D[TopSurf, r]

\(The\ oriented\ surface\ patch\ d \( 
          A\& -> \)\ for\ the\ top\ surface\ is\ the\ following\ \((multiplied\ by\ dr\ dθ)\)\)

NormalVecTopSurf  =  
      FullSimplify[Cross[TopSurfRr, TopSurfRt]]

TopIntegrandθr   =  
      FullSimplify[(FVecLondon /. {ξ -> 
                  r Cos[t],  η ->  r Sin[t], ζ -> 
                  L/2}) . NormalVecTopSurf]

Because we will need to numerically integrate over r and over t, we might try the same trick as above and see if one of these integrals can be expressed in closed form.

Integrating over θ doesn't provide a closed form (this next step takes a while to evaluate...)

Integrate[TopIntegrandζr,  {t, 0, 2  π}, 
      Assumptions  ->  
        r  ≥  0  &&  
          L  >  0  &&   x ∈  Reals  &&  
          y  ∈  Reals  &&  z  ∈  Reals]

Integrating over r does work

TopIntegrandθIndr =  
      Integrate[TopIntegrandθr,  r, 
        Assumptions  ->  
          t  ≥  0  &&  
            L  >  0  &&    x ∈  Reals  &&  
            y  ∈  Reals  &&  z  ∈  Reals]

TopIntegrandθ[x_, y_, z_, CylRad_, CylLen_]  := 
      Evaluate[Simplify[(TopIntegrandθIndr /. {r -> CylRad, 
                  L -> CylLen}) - (TopIntegrandθIndr /. {r -> 0, 
                  L -> CylLen})]]

TopContribution[xpos_,  ypos_,  zpos_,  CylRad_ ,  CylLen_]  :=  
      NIntegrate[
        Evaluate[
          TopIntegrandθ[xpos, ypos, zpos,  CylRad,  CylLen], {t, 0, 
            2 π}]]

TopContribution[dist_,  zpos_,  CylRad_ ,  CylLen_]  :=  
      NIntegrate[
        Evaluate[
          TopIntegrandθ[dist, 0, zpos,  CylRad,  CylLen], {t, 0, 
            2 π}]]

TopContribution[Sqrt[2], Sqrt[2], 1/3, 1, 2]

TopContribution[2, 1/3, 1, 2]

Bottom Cap Surface normals and differential quantities

BotSurf  =  {r Cos[t], r Sin[t], ((-L)/2)}

BotSurfRt  =  D[BotSurf, t]
BotSurfRr    =  D[BotSurf, r]

\(The\ oriented\ surface\ patch\ d \( 
          A\& -> \)\ for\ the\ bottom\ surface\ is\ the\ following\ \((multiplied\ by\ dr\ dθ)\)\)

NormalVecBotSurf  =  
      FullSimplify[Cross[BotSurfRt, BotSurfRr]]

BotIntegrandθr  =  
      FullSimplify[(FVecLondon /. {ξ -> 
                  r Cos[t],  η ->  
                  r Sin[t], ζ -> -L/2}) . 
          NormalVecBotSurf]

BotIntegrandθIndr =  
      Integrate[BotIntegrandθr,  r, 
        Assumptions  ->  
          t  ≥  0  &&  
            L  >  0  &&    x ∈  Reals  &&  
            y  ∈  Reals  &&  z  ∈  Reals]

BotIntegrandθ[x_, y_, z_, CylRad_, CylLen_]  := 
      Evaluate[Simplify[(BotIntegrandθIndr /. {r -> CylRad, 
                  L -> CylLen}) - (BotIntegrandθIndr /. {r -> 0, 
                  L -> CylLen})]]

BotContribution[xpos_,  ypos_,  zpos_,  CylRad_ ,  CylLen_]  :=  
      NIntegrate[
        Evaluate[
          BotIntegrandθ[xpos, ypos, zpos,  CylRad,  CylLen], {t, 0, 
            2 π}]]

BotContribution[dist_,  zpos_,  CylRad_ ,  CylLen_]  :=  
      NIntegrate[
        Evaluate[
          BotIntegrandθ[dist, 0, zpos,  CylRad,  CylLen], {t, 0, 
            2 π}]]

BotContribution[Sqrt[2], Sqrt[2], 1/3, 1, 2]

BotContribution[2, 1/3, 1, 2]

Adding up all the surface integral contributions

LondonCylinderPotential[dist_,  zpos_,  CylRad_ ,  CylLen_]   := 
      CylinderContribution[dist, zpos, CylRad, CylLen]  +  
        TopContribution[dist, zpos, CylRad, CylLen]  +  
        BotContribution[dist, zpos, CylRad, CylLen]

LondonCylinderPotential[2,  .5, 1, 3]

LondonCylinderPotential[ .25,  .5, 1, 3]

<< Graphics`Graphics`

Visualize result as a function of radial distance at different altitudes

LondonPlot  = 
      Plot[
{LondonCylinderPotential[dist, 0, 1, 4/3], 
          LondonCylinderPotential[dist, 2/3, 1, 4/3], 
          LondonCylinderPotential[dist, 4/3, 1, 4/3], 
          LondonCylinderPotential[dist, 2, 1, 
            4/3]}, 
{dist, 0.01, 3}, 
        PlotStyle -> {{Thickness[0.02], 
              RGBColor[1, 0, 0]}, {Thickness[0.015],  
              RGBColor[0, 0.5, 0]},  {Thickness[0.01], 
              RGBColor[0, 0, 1]}, {Thickness[0.005], 
              RGBColor[1, 0, 1]}}]

The plot above shows the London potential vs. radial distance for four different z values, z = 0, z = 2/3, z = 4/3, and z= 2. The cylinder height is 4/3 and its radius is 1; for values of z that intersect the cylinder, the potential diverges as dist -> 1.

Show[LondonPlot, PlotRange -> {-5, 3}]

The plot above shows the London potential vs radial distance for the same z values as in the previous plot, but for a much smaller range of potential values (-1, 0), to show the z dependence of the potential more clearly.

Visualize potential at fixed altitudes above cylinder

TopOfCylinder  = 
      Plot[{LondonCylinderPotential[dist, 1.1, 1, 1], 
          LondonCylinderPotential[dist, 1.2, 1, 1], 
          LondonCylinderPotential[dist, 1.3, 1, 1], 
          LondonCylinderPotential[dist, 1.4, 1, 1]}, {dist, 0, 3}, 
        PlotStyle -> {{Thickness[0.02], 
              RGBColor[1, 0, 0]}, {Thickness[0.015],  
              RGBColor[0, 0.5, 0]},  {Thickness[0.01], 
              RGBColor[0, 0, 1]}, {Thickness[0.005], 
              RGBColor[1, 0, 1]}}]

The plot above shows the London potential vs. radial distance for z values equal to and larger than the cylinder height. The cylinder height is 1 and its radius is 1;.the potential is plotted for z = 1.1, z = 1.2, z = 1.3, and z= 1.4.  All these z values are beyond the end of the cylinder, and the potentials do not diverge at dist = 1. Note that unlike the previous plots, this one extends from dist = 0 to 3.

The contour plot below would take an enormously long time to compute if we had not employed all of the ``integral tricks''

ContourPlot[
      LondonCylinderPotential[dist, height, 1, 0.25], {dist, 0.001, 
        2}, {height, 0.001, 2}, Contours -> 25, 
      ColorFunction -> (Hue[0.6 #] &)]