Generalizations of the Fundamental Theorem of Calculus (Part I)

Using Green's theorem in the plane to make the integration of the triangular charged surface patch simpler:

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)= \(TraditionalForm\`\(dξ\ dη\)\/\@\(\((x - ξ)\)\^\(\(2\)\(\ \)\) + \((y - η)\)\^\(\(2\)\(\ \)\) + z\^\(\(2\)\(\ \)\)\)\)

r[x_ ,  y_ ,  z_ ]  :=  
        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)\)\  = \(\(\@\(\(3\)\(\ \)\)\) t\)\/2)\)\ ; \ \((0\  < \ t\  < \ 1)\)\)
\(The\ other\ two\ sides\ \(are : \ \((ξ \((t)\)\  = \(-\(\(\(\ \)\(t\)\)\/2\)\)\ \ , \ η \((t)\)\  = \(\@\(\(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  (∫_▲)∫\(\((∂F\_2\/∂x\  - ∂F\_1\/∂y)\)\ dx dy\) = (∫_△) (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

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.

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

Doing the same integral numerically

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

ContourPlot[Pot[a, b, 1/10], {a, -1, 1}, {b, - .5, 1.5}, 
      Contours -> 8, 
      ColorFunction -> (Hue[1  - 0.66*#] &) ]

Surfaces

Now we look at several different ways in which surfaces can be represented algebraically...

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

GraphFunction[x_,  
        y_]  :=  (((x - y) (x + y))/(1  +  (((x + y))^2)))

Plot3D[
      Evaluate[GraphFunction[x, y]], {x, -1, 1}, {y, -1, 
        1}]

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

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

Here is a view of the same function from a different perspective...

ParametricPlot3D[
      Evaluate[SurfaceParametric[u, v]], {u, -2, 2}, {v, -2, 2}, 
      ViewPoint -> {1.632,  0.625,  -0.972}]

The next statement plots the same parametric function over a somewhat larger range of u and v, using a clever strategy that gradually builds up the surface and serves to aid visualization of itss complex shape. After computing the sequence of plots, select them and close the group, then animate the collection.

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

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

<< Graphics`ContourPlot3D`

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

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

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

cpa  = 
      ContourPlot3D[
        (x^2)  - 4 x y  +  (y^2)  +  y x (z^2), {x, -3, 
          3}, {y, -3, 3}, {z, -3, 3}, Contours -> {0, 2, 8}, 
        PlotPoints -> {5, 7}, 
        DisplayFunction -> Identity]

cpb  =  
      ContourPlot3D[
        (x^2)  - 4 x y  +  (y^2)  +  y x (z^2), {x, -3, 
          3}, {y, -3, 3}, {z, -3, 3}, Contours -> {0, 2, 8}, 
        PlotPoints  ->  {5, 7}, 
        ContourStyle ->  {{Hue[0]}, {Hue[ .25]}, {Hue[ .5]}}, 
        Lighting -> False, DisplayFunction -> Identity]

Show[GraphicsArray[{cpa, cpb}]]

Table[
      ContourPlot3D[
        (x^2)  - 4 x y  +  (y^2)  +  y x (z^2), {x, -3, 
          3}, {y, -3, 3}, {z, -3, 3}, Contours -> {i}, 
        PlotPoints -> {5, 7}], {i, -2, 10,  .5}]

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).)

spheresurf[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[spheresurf[u, v], u] // Simplify
Rv[u_ ,  v_]  =  
      D[spheresurf[u, v], v] // Simplify

<< Calculus`VectorAnalysis`

The normal vector:

NormalVector[u_ ,  v_ ]  =  
      CrossProduct[Ru[u, v], Rv[u, v]] // Simplify

NormalMag =  
      Sqrt[DotProduct[NormalVector[u, v], NormalVector[u, v]] // Simplify] // 
        PowerExpand

The Unit Normal: No surprise here, it is just the vector from the origin to a point on a unit sphere...

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

<< Graphics`ParametricPlot3D`

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, 
        2  Pi}, {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\ 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.