Generalizations of the Fundamental Theorem of Calculus (Part I)

Using Green's theorem to Reduce a Surface Integration to a Path Integral

Suppose there is a uniformly charged surface  (σ≡charge/area=1) occupying an equilaterial triangle in the z=0 plane:

The electrical potential goes like  (1/r), therefore the potential of a unit charge located at (x,y,z) from a small surface patch at (ξ,η,0) is ((σ dξ dη)/r)= ((dξ dη)/(Sqrt[(((x-ξ))^(2 ))+(((y-η))^(2 ))+(z^(2 ))]))

Try to convert the integral over the region to an integral over its boundary.  We may be able to find a closed form—and even if not the numerical integration should be much more efficient...

We will use Green's Theorem in the Plane:
∫∫(((∂(Subscript[F, 2]))/(∂ x))-((∂(Subscript[F, 1]))/(∂ y)))dxdy = ∫((Subscript[F, 1])dx + (Subscript[F, 2])dy)

Find a function, that when we take the derivative with respect to x, we get 1/r.
We will let (Subscript[F, 2]) be the function we find, and set (Subscript[F, 1])=0

The boundary defined by the upper right side of the triangle is (ξ(t) =(((1-t))/2), η(t) =(((Sqrt[3 ])t)/2)) ; (0 < t < 1)
The other two sides are: (ξ(t) =-(( t)/2)  , η(t) =(Sqrt[3 ])(((1-t))/2) )  (0 < t < 1) and
(ξ(t) =((-1)/2) + t  , η(t) = 0)   (0 < t < 1)

To use Green's theorem in the plane  (∫_▲)∫(((∂(Subscript[F, 2]))/(∂x)) -((∂(Subscript[F, 1]))/(∂y))) dxdy = (∫_△) (Subscript[F, 1])dx + (Subscript[F, 2])dy
Try to find a vector function ((Subscript[F, 1]), (Subscript[F, 2])) that matches the integrand over the domain.  Let (Subscript[F, 2]) = 0, then to find (Subscript[F, 1 ])which when differentiated by η gives 1/r

F1[x_,y_ , z_] =-Integrate[(1/(Sqrt[(((x-ξ))^2)+(((y-η))^2)+(z^2)])),η]

The third (horizontal) boundary of the triangle patch looks like the easiest, let's see if an integral can be found over that patch:

Bottomside = F1[x,y,z]/.{ξ->t-(1/2), η->0}//Simplify

The following attempt to find a closed-form integral seems to take "forever" to evaluate and I ended up aborting it, then continuing to evaluate the remainder of the notebook. This command has been commented out to avoid evaluation

(*
BottomContribution = Integrate[Bottomside,{t,0,1},Assumptions->x∈Reals&&y∈ Reals&&z>0]
*)

NEside = F1[x,y,z]/.{ξ->((1-t)/2), η->(((Sqrt[3]) t)/2)}//Simplify

NWside=F1[x,y,z]/.{ξ->((-t)/2), η->(((Sqrt[3]) (1-t))/2)}//Simplify

(∫_△) (Subscript[F, 1])dx =(∫_0^1) NEside ((-dt)/2) + (∫_0^1) NWside ((-dt)/2) + (∫_0^1) Bottomside dt

integrand =Simplify[ ((-(NEside + NWside))/2) +Bottomside]

The following attempt to find a closed-form integral seems to take "forever" to evaluate and I ended up aborting it, then continuing to evaluate the remainder of the notebook. This command has been commented out to avoid evaluation

(*
PotXYZ = Integrate[integrand,{t,0,1},Assumptions->x∈Reals&&y∈ Reals&&z>0]
*)

The Potential over a Triangular Charged Surface Patch: Example of Green's theorem and Faster Integration

Doing the same integral as in the previous lecture numerically, but this time over the boundary of the triangle instead of the triangle area.

Pot[X_, Y_, Z_] := NIntegrate[Evaluate[integrand/.{x->X,y->Y,z->Z}],{t,0,1}]

We will create contourplots (level sets of constant potential) at as a function of different heights. We check the timing of the computation to compare to method in the last lecture.

ncplot[h_] := ncplot[h] =ContourPlot[Pot[a,b,h],{a,-1,1},{b,-.5,1.5},Contours->Table[v,{v,.25,2,.25}],ColorFunction->ColorData["TemperatureMap"],ColorFunctionScaling->False,PlotPoints->11 ]
Timing[ncplot[.05]]

Row[{TextCell["Computing ContourPlots a different h: Progress:  ","Text"],ProgressIndicator[Dynamic[h],{0,.5}]}]
ncplots = Table[ncplot[h],{h,.025,.5,.025}];

ListAnimate[ncplots]

Representations of Surfaces (part 1):  z=f(x,y)

Surface of the form: (z=f(x,y))

We construction an example function of the form z = U(S,V) for display with different kinds of superimposed information.

GraphFunction[x_, y_] := (((x-y)(x+y))/(1 + (((x+y))^2)))
assump = {x ∈ Reals, y∈Reals};

A default view of the surface.

plotdefault=Plot3D[GraphFunction[x,y],{x,-3,3},{y,-3,3},PlotLabel->"Default"]

Plotting the surface with levels of constant height superimposed and colored according to its height.

plotlevels=Plot3D[GraphFunction[x,y],{x,-3,3},{y,-3,3},MeshFunctions->( #3&),ColorFunction->"Rainbow",PlotLabel->"Constant Heights"]

This function will be used to paint surface according to an angular variable.

angle[x_] := ((Pi/2 + ArcTan[x])/Pi)
angle[x_, y_] := ((Pi/2 + ArcTan[x,y])/Pi)

Drawing circles on the plot and coloring accoding to the angle; this is superposing a visualization of the cylindrical coordinate system.

plotcircles =Plot3D[GraphFunction[x,y],{x,-3,3},{y,-3,3},MeshFunctions->(Sqrt[#1^2 + #2^2]&),ColorFunction->(Hue[angle[#1,#2]*0.5]&),ColorFunctionScaling-> False,PlotLabel->"Cylindrical Coordinates"]

Here we write a function to calculate the curvature of f(x,y) as a function of x and y. We will use this to draw lines of constant curvature.

CurvatureOfGraph[f_, x_, y_] := FullSimplify[Module[{dfdx=D[f[x,y],x],dfdy=D[f[x,y],y],d2fdx2 = D[f[x,y],{x,2}],d2fdy2 = D[f[x,y],{y,2}],d2fdxdy= D[f[x,y],x,y]} ,Return[((1 + dfdx^2) d2fdx2 - 2 dfdx dfdy d2fdxdy + (1+ dfdy^2)d2fdy2)/Sqrt[1 + dfdx^2 + dfdy^2]]],Assumptions->assump]

This is a moderately advanced application---here I define a function symbolically by evaluating the curvature function, this will be passed to Plot3D as the MeshFunction.

CurvFunc = Function[{x,y},Evaluate[CurvatureOfGraph[GraphFunction,x,y]]]

Representations of Surfaces (part 2):  z=f(x,y)

dfdx = Function[{x,y},Evaluate[FullSimplify[D[GraphFunction[x,y],x],Assumptions->assump]]]
dfdy = Function[{x,y},Evaluate[FullSimplify[D[GraphFunction[x,y],y],Assumptions->assump]]]

This is the surface with lines of constant curvature superimposed, and with colors associated with the local normal.

plotcurvature =Plot3D[GraphFunction[x,y],{x,-3,3},{y,3,-3},MeshFunctions->(CurvFunc[#1,#2]&),
MeshStyle->Thick,PlotLabel->"Curvatures(level sets) and Normals(color variation)",ColorFunction->(Glow[RGBColor[angle[dfdx[#1,#2]],angle[dfdy[#1,#2]],0.75]]&),ColorFunctionScaling->False,Lighting->None]

Visualizing all the examples together.

GraphicsGrid[{{plotdefault,plotlevels},{plotcircles,plotcurvature}},ImageSize->Full]

Representations of Surfaces:  z=f(x,y) (Frivolous—and somewhat advanced—Examples) Floating Pixels

MinMax[alist_List] := Module[{flatlist=Flatten[alist]},Return[{Min[flatlist],Max[flatlist]}]]
mug =Import["http://pruffle.mit.edu/~ccarter/ch_face_frames/Carter_2000_verysmall.png"];
ProgressIndicator[Dynamic[i],{1,64}]
vp[i_] := {.1 Sin[(i-1)  Pi/31], Sin[(i-1) 2 Pi/31],2Cos[2(i-1) Pi/63]};
minmax=MinMax[mug[[1,1]]];
Table[mugshot[i]=ListPlot3D[mug[[1,1]],MeshStyle->None, Mesh->None,InterpolationOrder->0,ColorFunction->"GreenBrownTerrain",Axes->False,ViewPoint->vp[i],
PlotRange->minmax,
ImageSize->Full,SphericalRegion->True];,{i,1,64}];
Manipulate[mugshot[frame],{frame,1,64,1}]

Representations of Surfaces:  z=f(x,y) (Frivolous—and somewhat advanced—Examples) Morphing and Warping Images

Grid[{{TextCell["Reading Files off Net","Text"],ProgressIndicator[Dynamic[read],{1,28}]},
{TextCell[ "Graphics Computation","Text"],ProgressIndicator[Dynamic[gcomp],{1,28}]}}]
Table[chface[read] =Import["http://pruffle.mit.edu/~ccarter/ch_face_frames/ch_face."<>ToString[100 + read-1]<>".png"];facedata[read] =ListInterpolation[chface[read][[1,1]],{{0,1},{0,1}}];
If[read==1,minmax= MinMax[chface[read][[1,1]]];,minmax =  MinMax[{minmax,chface[read][[1,1]]}]];,{read,1,28,1}];
pface[i_] :=  Plot3D[facedata[i][x,y],{y,0,1},{x,0,1},PlotRange->minmax,ColorFunction->"GreenBrownTerrain",Mesh->False,Axes->False,ViewPoint->{-0.25,-2,5},ImageSize->Full]
ListAnimate[Table[pface[gcomp], {gcomp,1,28,1}],DefaultDuration->10]

Surface of the form: (x(u,v), y(u,v), z(u,v))

SurfaceParametric[u_, v_] := {Cos[u] v, u Cos[u+v], Cos[u]/(.1 + Cos[u]^2)}

ParametricPlot3D[Evaluate[SurfaceParametric[u,v]],{u,-2,2},{v,-2,2}]

Using Manipulate, we can vary the boundary domain, and provide a more intuitive way to understand this complicated surface.

evolution=Table[ParametricPlot3D[Evaluate[SurfaceParametric[u,v]],{u,-ep,ep},{v,-ep,ep},PlotRange->{{-4,4},{-4,4},{-4,4}},PlotPoints->{1+Round[ep/.125],1+Round[ep/.125]},ImageSize->Full],{ep,.125,4.25,.125}];
ListAnimate[evolution,ImageSize->Full]

Surface of the form: (F(x,y,z) = constant)

ConstFunction =  (x^2) -4 x y + (y^2) + (z^2)

ContourPlot3D[ConstFunction,{x,-1,1},{y,-1,1},{z,-1,1},Contours->{2.5}]

The following statements produce contour plots of the same function, using two different methods for colorizing the surfaces...

cpa =ContourPlot3D[ConstFunction,{x,-3,3},{y,-3,3},{z,-3,3},Contours->{0,2,8}]

cpb =ContourPlot3D[ConstFunction,{x,-3,3},{y,-3,3},{z,-3,3},Contours->{0,2,8},ContourStyle->{
Directive[Pink,Opacity[0.8]],Directive[Yellow,Opacity[0.8]],Directive[Orange,Opacity[0.8]]}]

Manipulate[ContourPlot3D[ConstFunction,{x,-3,3},{y,-3,3},{z,-3,3},Contours->{i},ImageSize->Full],{i,-2,10,.25}]

Surface Integration

Example: Integrating an Orientation-Dependent Surface Tension over the surface of a cube and a sphere

Vector function for surface of sphere of radius R, u is longitude (u ∈ (0, 2π)) and v is latitude (v ∈ (-(π/2), (π/2)))
(Note: this is the parametric form of the equation of the sphere (x^2)+(y^2)+(z^2)=(R^2).)

sphere[u_ , v_ ] := R{Cos[v] Cos[u] , Cos[v] Sin[u] , Sin[v]}

The coordinate vectors at (u,v) (Note: these vectors are tangent to the sphere surface at the point (u,v).)

Ru[u_ , v_] = D[sphere[u,v],u]//Simplify
Rv[u_ , v_] = D[sphere[u,v],v]//Simplify

A normal vector and its unit normal: no surprise here, the unit normal is just the vector from the origin to a point on a unit sphere...

Needs["VectorAnalysis`"]
NormalVector[u_ , v_ ] = CrossProduct[Ru[u,v],Rv[u,v]]//Simplify
NormalMag=FullSimplify[ Norm[NormalVector[u,v]],Assumptions->{R ≥ 0, 0 ≤ u ≤ 2 π,-π/2 < v < π/2}]
UnitNormal[u_, v_] = NormalVector[u,v]/NormalMag

Now, suppose the surface energy ( (Joules/(meter^2))) depends on the local orientation of the surface. This is case for many single crystals. If the single crystal is cubic, it is often the case that the (100) planes have lower energy than the (111) planes.  The following model for orientation-dependent surface tension has this behavior and can be employed for illustrative purposes:

SurfaceTension[nvec_] :=1 + (Subscript[gamma, 111])*((nvec[[1]])^2)((nvec[[2]])^2) ((nvec[[3]])^2)

Here is a plot of our model aniostropic surface energy-- the lobes at [111] indicate that those orientations have higher energy:

SphericalPlot3D[SurfaceTension[UnitNormal[u,v]]/.(Subscript[gamma, 111]) -> 12,{u,0,2Pi},{v,-Pi/2,Pi/2}]

The total energy of a spherical particle of unit radius having the anisotropic surface tension represented above is given by the double integral

SphereEnergy =Integrate[Integrate[SurfaceTension[UnitNormal[u,v]] Cos[v],{u,0,2π}],{v,-(π/2),(π/2)}]

Compare to Cube of same volume as the unit sphere: (CubeSide^3) = ((4 π)/3)

CubeSide = (((4 π)/3)^(1/3))

CubeEnergy = 6 ( (CubeSide^2) SurfaceTension[{1,0,0}])

Find the condition on (Subscript[gamma, 111]) such that the two surfaces have the same energy

EqualEnergies = Solve[CubeEnergy==SphereEnergy, (Subscript[gamma, 111])]//Flatten

If the anisotropy factor is greater than about 25, then a cube is much lower energy

N[(Subscript[gamma, 111])/.EqualEnergies]

This computation compared two fixed shapes of the same volume.  To find which shape of all the possible shapes, one would use a procedure known as the Wulff construction.