Generalizations of the Fundamental Theorem of Calculus (Part II)
(note: some of the symbolic integrations in this notebook will take a long time to evaluate)

Using Divergence theorem to find Hamaker Field (here modeled as point-wise(1/(r^6))) for a Finite Body

The induced dipole—dipole interaction goes as -1/(r^6), and it is this term which shows up in the Lennard-Jones potential.  This dipole-dipole interaction gives rise to the Hamaker forces between two bodies.  As a first approximation to the interaction between two bodies, one must intergrate 1/(r^6) from a single point in the first body to all points in the second body.  The net interaction energy is the sum over all such points and therefore this energy produces a six-dimensional integration. There are a few cases, such a two spheres, a sphere and a seminfinite half-space, that can be worked out exactly.  However, to model the interaction energy of more complex shapes, there are no closed-form solutions and one must resort to numerical integration.

Numerical integration is a cpu-time-consuming numerical procedure—a six dimensional integration can be very expensive (time is money).  If there is a way to reduce the dimensionality of the integration, then we can reap rewards for our cleverness.  One trick is to use the divergence theorem to push the integration over a volume, to an integration over a surface. For example, we could use the divergence theorem:
∫∫∫(Subscript[ , volume]) ∇· (OverVector[P]) dV = ∫∫(Subscript[ , surface])(OverVector[P])·d(OverVector[A])
We will work out such an example for a cylindrical volume interacting with a point in space, in other words, the London Interaction Potential for a finite cylinder. We will perform as many exact integrals as possible and reduce the cost of the numerical integration.  To do this, we must find a vector potential (OverVector[P]) such that ∇· (OverVector[P]) = -1/|(OverVector[r])-(OverVector[x])((|)^6) where (OverVector[r]) is a position in the integrated volume and  (OverVector[x]) is a point at which the potential is measured.  (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.  We use (CX,CY,CZ) to represent points in the cylinder and (X,Y,Z) to be points in space

PVecLondon = (1/(3((((((CX-X))^2) +(((CY-Y))^2) + (((CZ-Z))^2)))^3))){CX-X, CY-Y, CZ-Z}

Needs["VectorAnalysis`"]

The following verifies that the correct vector potential is obtained, the divergence must be taked with respect to the coodinates of the integration:

FullSimplify[Div[PVecLondon,Cartesian[CX,CY,CZ]]]

We will integrate over a cylinder of radius R and length L along the z-axis, with its middle at the origin.  First, let's use the radius of the cylinder to scale all the length variables: Let (X,Y,Z)/R = (x,y,z); (CX,CY,CZ)/R = (cx,cy,cz), and L/R = λ (the cylinder's aspect ratio).

ScaleRules={X->x R, Y-> y R, Z-> z R, CX-> cx R, CY-> cy R, CZ-> cz R};
PvecR5 = FullSimplify[R^5 PVecLondon/.ScaleRules,Assumptions->R> 0]

Therefore, ϕ((OverVector[x])) =  ∫∫∫(Subscript[ , volume]) ((-1)/((((Overscript[r, ->])-(Overscript[x, ->])))^6)) dV = ∫∫(Subscript[ , surface])(OverVector[PVecLondon])·d(OverVector[A])  = (∫∫(Subscript[ , cylinder
surface])   (OverVector[PVecR5])•d(Overscript[A, ->]) + ∫∫(Subscript[ , cylinder
ends])     (OverVector[PVecR5])•d(Overscript[A, ->]))/(R^5)
 is the total interaction between a point an a cylinder.  We can exploit the symmetry of the cylinder: we know that the potential will be indendent of the angular coordinate in cylindrical coordinates; it will be a function of ρ = (Sqrt[x^2+ y^2])and z.

We will do three integrals over the cylindrical surfaces using this expression to define the cylinder: (cx,cy,cz) = (Cos[θ], Sin[θ], cz):
The cylindrical surface  is the domain  θ ∈ (0, 2π), cz ∈ (-(λ/2) , (λ/2))
The two caps r ∈ (0,1),  θ ∈ (0, 2π), cz=±(λ/2)

Cylinder Surface Normals and Integrands

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 = { Cos[θ], Sin[θ],cz}

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

CylSurfRθ = D[CylSurf,θ]
CylSurfRcz  = D[CylSurf,cz]

The surface normal given by (OverVector[(Subscript[R, u])]) × (OverVector[(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[CylSurfRθ,CylSurfRcz]

The integrand to be evaluated over the cylinder surface is the vector potential, dotted into the normal vector.  Because of the cylindrical symmetry of this model, we can convert to cylindrical coordinates.  One set of coordinates is for the cylinder surface (ξ -> R Cos[θ], η -> R Sin[θ]) for fixed radius R (which is a model parameter) and another set of coordinates for where we will be testing the potential (x -> ρ Cos[α], y -> ρ Sin[α]).  Because the potential must be independent of α, we might as well set it to zero.

CylinderIntegranddθdζ = FullSimplify[(PvecR5/.{cx->Cos[θ], cy-> Sin[θ],x-> ρ ,y->0}).NormalVecCylSurf]

We have a choice whether to integrate over cz ∈  (-(λ/2) , (λ/2))  or  θ ∈ (0, 2π) first.  If we can a closed form for the cylinder surface over cz and then the cylinder end over r, then we can integrate the sum of these over θ together.

Try and see if we can do one of the two integrals---we have a choice of integrating over θ or (ζ for the cylinder sides, and  R) for the cylinder ends. We find a closed form for integrating  ζ for the sides and R for the top, and then subsequently numerically integrate θ for (0,2 π).


This will take a while to finish, even on a very fast  machine (as of 2007, one minute to do the integral and about 6 hours minutes to simplify the expression).
For this reason, I' ve set up an option to read the result over the web, instead of having each of your computers spend time on a previously computed result:
(NOTE: As of Mathematica 6.01 Oct. 2007,  If we first try the definite integral over θ ∈ (0, 2 π) , we get  ZERO which is  an incorrect result (I will include a discussion of why it is incorrect in the future))

Integrating over the Cylinder's Surface

UpperPlane={λ> 0,ρ>0, z>0 ,0 < θ< 2 π};

CylinderIntegrandUpperZdθ=
If[ChoiceDialog["Evaluate on your Computer (6 hours on a 2×2.66 GHz Dual-Core Intel Xeon) or Read a Precomputed Result From URL (a few seconds if you have a net connection)?",
{"Read URL"->True,"Evaluate Yourself"->False}],Import["http://pruffle.mit.edu/3.016-2007/CylinderIntegrandUpperZdtheta.m"],FullSimplify[Integrate[CylinderIntegranddθdζ,{cz,-λ/2,λ/2} , Assumptions->UpperPlane],Assumptions->UpperPlane]]

Here we restrict z to the upper half-space. We will treat z=0 below.

Here is the limit of  the integral for z> 0 (CylinderIntegranddθ) in the limit as z  -> 0.

CylinderIntegranddθZeroLimit =FullSimplify[Limit[CylinderIntegrandUpperZdθ,z-> 0]]

Here is the limit of  the integral z=0, it is not obvious that the limit and its value at z=0 are the same. (the following evaluation will take about 10-20 seconds)

FewerAssumptions= { R > 0 ,  λ> 0,ρ> 0, 0< θ < 2 π};

CylinderIntegrandAtZerodθ=Integrate[Evaluate[CylinderIntegranddθdζ/.z->0],{cz,-λ/2,λ/2} , Assumptions->FewerAssumptions]

The limit as z -> 0 and the integrand at z = 0 are the same, so we can use a single integrand

CylinderIntegranddθ[dist_,height_, AspectRat_]:= Evaluate[CylinderIntegrandUpperZdθ/.{ρ->dist,  z -> height,λ->  AspectRat}]
?CylinderIntegranddθ

CylinderContribution[dist_, height_, AspectRat_] := NIntegrate[CylinderIntegranddθ[dist,height, AspectRat],{θ,0,2 π}]

Integrating over the Cylinder Top

Now, perform the integral for the top and bottom surfaces

TopSurf = {r Cos[θ],r Sin[θ],(λ/2)}

TopSurfRθ = D[TopSurf,θ]
TopSurfRr   = D[TopSurf,r]

Now weThe oriented surface patch d(Overscript[A, ->]) for the top surface is the following (multiplied by dr dθ)

NormalVecTopSurf = FullSimplify[Cross[TopSurfRr,TopSurfRθ]]

EndAssumptions={λ> 0,ρ>0  ,z > 0, -1 ≤ Cos[θ]<1};
TopIntegranddθdr = FullSimplify[(PvecR5/.{cx->r Cos[θ], cy-> r Sin[θ],cz-> λ/2,x-> ρ,y->0}).NormalVecTopSurf,Assumptions->EndAssumptions] (* 5 seconds *)

If[ChoiceDialog["Evaluate on your Computer (6 minutes on a 2×2.66 GHz Dual-Core Intel Xeon) or Read a Precomputed Result From URL (a few seconds if you have a net connection)?",
{"Read URL"->True,"Evaluate Yourself"->False}],Import["http://pruffle.mit.edu/3.016-2007/TopIntegrands.m"],
TopIntegrandInsideAbovedr =Integrate[TopIntegranddθdr,{θ,0, 2 π}, Assumptions->{ 0 < r < 1,λ> 0,ρ<1  ,z > λ/2} ];
TopIntegrandInsideBelowdr =Integrate[TopIntegranddθdr,{θ,0, 2 π}, Assumptions->{ 0 < r < 1,λ> 0,ρ<1  ,z < λ/2} ];
TopIntegrandOutsideAbovedr=Integrate[TopIntegranddθdr,{θ,0, 2 π}, Assumptions->{ 0 < r < 1,λ> 0,ρ>1  ,z > λ/2} ];
TopIntegrandOutsideBelowdr=Integrate[TopIntegranddθdr,{θ,0, 2 π}, Assumptions->{ 0 < r < 1,λ> 0,ρ>1  ,z < λ/2} ]];
Grid[{{TopIntegrandInsideAbovedr,
TopIntegrandInsideBelowdr},
{TopIntegrandOutsideAbovedr,
TopIntegrandOutsideBelowdr}}]

TopIntegranddr=TopIntegrandInsideAbovedr

TopPart=If[ChoiceDialog["Evaluate on your Computer (6 minutes on a 2×2.66 GHz Dual-Core Intel Xeon) or Read a Precomputed Result From URL (a few seconds if you have a net connection)?",
{"Read URL"->True,"Evaluate Yourself"->False}],Import["http://pruffle.mit.edu/3.016-2007/TopPart.m"],Integrate[TopIntegranddr,{r,0,1},Assumptions->{ λ> 0,ρ>0  ,z > 0, z ≠λ/2}]]

TopContribution[dist_,height_, AspectRat_]:= Evaluate[TopPart/.{ρ->dist,  z -> height,λ->  AspectRat}]
?TopContribution

Integrating over the Cylinder Bottom

BotSurf = {r Cos[θ],r Sin[θ],((-λ)/2)}

BotSurfRθ = D[BotSurf,θ]
BotSurfRr   = D[BotSurf,r]

The oriented surface patch d(Overscript[A, ->]) for the bottom surface is the following (multiplied by dr dθ)

NormalVecBotSurf = FullSimplify[Cross[BotSurfRθ,BotSurfRr]]

BotIntegranddθdr = FullSimplify[(PvecR5/.{cx->r Cos[θ], cy-> r Sin[θ],cz-> -λ/2,x-> ρ,y->0}).NormalVecBotSurf,Assumptions->EndAssumptions]

If[ChoiceDialog["Evaluate on your Computer (2 minutes on a 2×2.66 GHz Dual-Core Intel Xeon) or Read a Precomputed Result From URL (a few seconds if you have a net connection)?",
{"Read URL"->True,"Evaluate Yourself"->False}],Import["http://pruffle.mit.edu/3.016-2007/BotIntegrandsdr.m"],
BotIntegrandInsidedr=Simplify[Integrate[BotIntegranddθdr,{θ,0, 2 π}, Assumptions->{ 0 < r < 1,λ> 0,ρ<1  ,z > 0} ],Assumptions->{ 0 < r < 1,λ> 0,ρ<1  ,z > 0}];
BotIntegrandOutsidedr=Simplify[Integrate[BotIntegranddθdr,{θ,0, 2 π}, Assumptions->{ 0 < r < 1,λ> 0,ρ>1  ,z >0} ],Assumptions->{ 0 < r < 1,λ> 0,ρ>1  ,z >0} ];]
Column[{BotIntegrandInsidedr,BotIntegrandOutsidedr}]

BotIntegranddr= BotIntegrandOutsidedr

Next evaluation will take about 40 seconds

BotPart= Integrate[BotIntegranddr,{r,0,1},Assumptions->{ λ> 0,ρ>0  ,z > 0}]

BotContribution[dist_,height_, AspectRat_]:= Evaluate[BotPart/.{ρ->dist,  z -> height,λ->  AspectRat}]

LondonCylinderPotential[dist_, height_, AspectRat_]  :=CylinderContribution[dist,height,AspectRat] + TopContribution[dist,height,AspectRat] + BotContribution[dist,height,AspectRat]

Speed and Efficiency: When to Evaluate the Right-Hand-Side during a Function Definition

Let's set a baseline to check efficiency. Here we check timing to integrate something

Timing[Integrate[Exp[Tan[x]],{x,0,c}]]

We check the same thing again, because Mathematica may have spent some time loading algorithms to integrate

Timing[Integrate[Exp[Tan[x]],{x,0,c}]]

Here, we time how long it takes to create a function (with delayed assignment), but using Evaluate on the rhs.

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

The following is equivalent to the above (safer) definition---and will work so long as c is not assigned to an expression.

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

The following should take the *least* amount of time to perform, but as we shall see is not as efficient in the long run.

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

?DelayedEvaluated
?Immediate
?FunctionDef

The following should give a rapid result

Timing[DelayedEvaluated[0.5]]
Timing[Immediate[0.5]]

The following will not be rapid, because it has to do the symbolic integration before returning the result.

Timing[FunctionDef[0.5]]

Visualizing the London Dispersion Potential in the Vicinity of Finite Cylinder: ContourPlots

Plot[LondonCylinderPotential[x,0.01,1],{x,0.0,1.5}]

Plot[LondonCylinderPotential[x,0.01,1],{x,0.0,1.5},Exclusions-> {1.0}, PlotRange->{-50,20},PlotStyle->{Thick,Red},BaseStyle->{Medium},ImageSize->Full,PlotLabel->"Hamaker Potential for Cylinder with Aspect Ratio 4\n at midplane, Height=0",AxesLabel-> {"\!\(\*FractionBox[\(distance\\ from\\ cylinder\), \(cylinder\\ radius\)]\)","Potential (cylinder radius)^5"}]

Plot[LondonCylinderPotential[0.0,h,1],{h,0.0,1.5},Exclusions-> {0.5}, PlotRange->{-50,20},PlotStyle->{Thick,Darker[Blue]},BaseStyle->{Medium},ImageSize->Full,PlotLabel->"Hamaker Potential for Cylinder with Aspect Ratio 4\n along axis of cylinder",AxesLabel-> {"\!\(\*FractionBox[\(height\\ above\\ midplane\), \(cylinder\\ radius\)]\)","Potential (cylinder radius)^5"}]

Plot[Log[Abs[LondonCylinderPotential[0.0,h,4]]],{h,0,2.5},Exclusions-> {2.0}, PlotRange->{0,18},PlotStyle->{Thick,Darker[Blue]},BaseStyle->{Medium},ImageSize->Large,PlotLabel->"Hamaker Potential for Cylinder with Aspect Ratio 4\n along axis of cylinder",AxesLabel-> {"\!\(\*FractionBox[\(height\\ above\\ midplane\), \(cylinder\\ radius\)]\)","Potential (cylinder radius)^5"}]

Plot[Log[Abs[LondonCylinderPotential[x,0,1]]],{x,0.0,1.5},Exclusions-> {1.0}, PlotRange->{0,20},PlotStyle->{Thick,Red},BaseStyle->{Medium},ImageSize->Large,PlotLabel->"Hamaker Potential for Cylinder with Aspect Ratio 4\n at midplane, Height=0",AxesLabel-> {"\!\(\*FractionBox[\(distance\\ from\\ cylinder\), \(cylinder\\ radius\)]\)","Potential (cylinder radius)^5"}]

If[ChoiceDialog["Evaluate on your Computer (5 minutes on a 2×2.66 GHz Dual-Core Intel Xeon) or Read a Precomputed Result From URL (a few seconds if you have a net connection)?",
{"Read URL"->True,"Evaluate Yourself"->False}],Import["http://pruffle.mit.edu/3.016-2007/conplotouter.m"],
conplotouter = ContourPlot[-Log[Abs[LondonCylinderPotential[dist,h,4]]],{dist,0,1.5},{h,0,2.5},RegionFunction->Function[{dist,h},dist> 1.01 || h> 2.01],PlotRange->{-15,1},Contours->15, ColorFunction->"AvocadoColors",AspectRatio->Automatic,ImageSize->Medium,PlotLabel->"Hamaker Potential for \n Cylinder with Aspect Ratio 4\n-[Log[Abs[Potential (Radius^5)]] for exterior\nLog[Potential (Radius^5)] for interior]",BaseStyle->{Medium},AxesLabel-> {"(distance from cylinder)/radius","(height above midplane)/radius"},AspectRatio->Automatic,Exclusions-> {{dist==1.0,Abs[dist-1.0] < 0.001},{h==2.0, Abs[h-2] < 0.001}}];];
conplotouter

If[ChoiceDialog["Evaluate on your Computer (16 minutes on a 2×2.66 GHz Dual-Core Intel Xeon) or Read a Precomputed Result From URL (a few seconds if you have a net connection)?",
{"Read URL"->True,"Evaluate Yourself"->False}],Import["http://pruffle.mit.edu/3.016-2007/conplotinner.m"],
conplotinner = ContourPlot[Log[LondonCylinderPotential[dist,h,4]],{dist,0,1.5},{h,0,2.5},RegionFunction->Function[{dist,h},dist< 0.99 && h<1.99],PlotRange->{0,15},Contours->16, ColorFunction->"LakeColors",AspectRatio->Automatic,ImageSize->Medium,PlotLabel->"Hamaker Potential for \n Cylinder with Aspect Ratio 4\n-[Log[Abs[Potential (Radius^5)]] for exterior\nLog[Potential (Radius^5)] for interior]",BaseStyle->{Medium},AxesLabel-> {"(distance from cylinder)/radius","(height above midplane)/radius"},AspectRatio->Automatic];];
conplotinner

Show[conplotinner,conplotouter]

Visualizing the London Dispersion Potential in the Vicinity of Finite Cylinder: Three-Dimensional Graphics

plotoutside3D =
If[ChoiceDialog["Evaluate on your Computer (1 minute on a 2×2.66 GHz Dual-Core Intel Xeon) or Read a Precomputed Result From URL (a few seconds if you have a net connection)?",
{"Read URL"->True,"Evaluate Yourself"->False}],Import["http://pruffle.mit.edu/3.016-2007/plotoutside3D.m"],
Plot3D[-Log[Abs[LondonCylinderPotential[dist,h,4]]],{dist,0,1.5},{h,0,2.5},RegionFunction->Function[{dist,h},dist> 1.01 || h> 2.01],PlotRange->{-15,1},MeshFunctions->{#3&},ColorFunction->"AvocadoColors",AspectRatio->Automatic,ImageSize->Large,BaseStyle->{Medium},AxesLabel-> {"distance/radius","height/radius"},AspectRatio->Automatic]];
plotoutside3D

plotinside3D =
If[ChoiceDialog["Evaluate on your Computer (1 minute on a 2×2.66 GHz Dual-Core Intel Xeon) or Read a Precomputed Result From URL (a few seconds if you have a net connection)?",
{"Read URL"->True,"Evaluate Yourself"->False}],Import["http://pruffle.mit.edu/3.016-2007/plotinside3D.m"],
Plot3D[Log[Abs[LondonCylinderPotential[dist,h,4]]],{dist,0,1.5},{h,0,2.5},RegionFunction->Function[{dist,h},dist< 0.99 && h< 1.99],PlotRange->{0,15},MeshFunctions->{#3&},ColorFunction->"LakeColors",AspectRatio->Automatic,ImageSize->Large,BaseStyle->{Medium},AxesLabel-> {"distance/radius","height/radius"},AspectRatio->Automatic]];
plotinside3D

Show[plotinside3D,plotoutside3D,PlotRange->{-15,15}]