4020.sagews

Calculate the surface integral of flux type  $\iint_{S^+} P(x,y,z)\, dydz+Q(x,y,z) \,dxdz+R(x,y,z) \,dxdy$ where $S^+$ is a surface with given orientation.

This is equivalent to the problem: Find the flux of the vector field $\vec F(x,y,z)=\,\lt P,Q,R \gt\,$   over $S^+$.

We will solve the (relatively) simple problem of $\iint_S x dydz+y \,dxdz+z \,dxdy$ where S is top side of the triangle x+y+z=1, $x, \,y, \,z \ge 0$. Other examples follow.

---

Formula Sheet From the formula sheet, we see:  $\iint_S^+ P dydz+Q \,dxdz+R \,dxdy = \iint_S \vec F \cdot d \vec S$, where $\vec F= \lt P,Q,R \gt$

Surface element is the vector product:  $d \vec S= \vec n \cdot dS = \pm \,\, d \vec S_u \times d \vec S_v \,\,\, du \, dv$

Remember that the $\pm$ depends on whether or not the parametrization is orientation preserving. We will find this by finding the sign of the dot product of the vector: $\vec n$ and the vector: $d \vec S_u \times d \vec S_v$.

Integrand is the mixed product: $\vec F \,d \vec S = \pm \,\, \vec F \cdot d \vec S_u \times d \vec S_v \,\,\, du\, dv$

Preparation for SOLVER: We must parameterize the surface S with 2 parameters, find the intervals and find any vector $\vec n$ normal to S and in the positive direction of the given orientation.

Work particular to the given problem: 

== We have an explicit function $z=f(x,y)$ for the surface. Namely $z=1-x-y$ So we just let $x=u$, $y=v$ and $z=f(u,v)=1-u-v$. Other parametrizations will not be so easy.

== So: $\left\{ \begin{array}{l}x = u \\y = v\\z =1-u-v \end{array} \right.\,\,\,\,\,u \in [0,1 ] \,\,\,\,v \in [0,-u+1]$
   (To get the intervals, we look at the vertices of the triangle in the u0v plane. They are V1(u=1,v=0), V2(u=0,v=1) and V3(u=0,v=0). We draw this region. It is a triangle where v goes from 0 to the line joining V1 and V2, namely v=-u+1 and then u goes from 0 to 1.)

== By examination, we see: $\vec n= <1,1,1>$   (That is, the vector pointing up from the triangle.)

---

SOLVER: Our program follows the hand solution method. Everything in red requires something particular to your problem!

1. Parameterize S and input as vector function. Requires parametrization - done above.
2. Find $\vec n$ and input as vector function. Requires normal - done above.
3. Input F.
4. Find the partials of S. 
5. Determine $\pm$ by checking orientation preserving (+1 if orientation-preserving, else -1). Requires checking that you get $\pm 1$.
6. Substitute parameterization in F.
7. Find the mixed product  $\vec F \,d \vec S = \pm \,\,\,\vec F \cdot d \vec S_u \times d \vec S_v \,\,\, du\, dv$
8. Find intervals of integration and integrate. Requires intervals of parametrization - done above.

---

Step 0: The program declares our variables. We are given variables (x,y,z). We need (u,v) as parameters.

var ('u v'); var('x y z')

\newcommand{\Bold}[1]{\mathbf{#1}}\left(u, v\right)

\newcommand{\Bold}[1]{\mathbf{#1}}\left(x, y, z\right)

---

Step 1: We define $\vec S$, $\vec n$ and $\vec F$. Changes with problem; you must input your parametrization in Sf and your vector in n.

Sf=vector((u,v,1-u-v))
n=vector((1,1,1))
F=vector((x,y,z))

---

Step 2: The program finds the partial derivatives.

Sprime_u=diff(Sf,u)
view(Sprime_u)

\newcommand{\Bold}[1]{\mathbf{#1}}\left(1,\,0,\,-1\right)

Sprime_v=diff(Sf,v)
view(Sprime_v)

\newcommand{\Bold}[1]{\mathbf{#1}}\left(0,\,1,\,-1\right)

---

Step 3: The program checks whether our parameterization is "orientation-preserving" or not. Check that we get 1 or -1.

If we get 0 here, we need to change the point, e.g. (u=2.0,v=2.0). Use "real" values with decimal points. (If we get something other than +1, -1 or 0, we have made an error in step 1.)

npar=Sprime_u.cross_product(Sprime_v)
np=sign((n.dot_product(npar))(u=1.0,v=1.0))
view(np)

\newcommand{\Bold}[1]{\mathbf{#1}}1

---

Step 4: The program defines a standard sage function for "change of variables" for a 2-variable parametrization (surface).

The program changes the variables of $\vec F$.

def changevar(f, eqn, newvar1,newvar2):
    return f.substitute(eqn)

F=changevar(F,x==Sf[0],u,v)
F=changevar(F,y==Sf[1],u,v)
F=changevar(F,z==Sf[2],u,v)
view(F)

\newcommand{\Bold}[1]{\mathbf{#1}}\left(u,\,v,\,-u - v + 1\right)

---

Step 5: Then program finds the mixed product and multiplies it by the orientation from step 3.

M=matrix(([F[0],F[1],F[2]],[Sprime_u[0],Sprime_u[1],Sprime_u[2]],[Sprime_v[0],Sprime_v[1],Sprime_v[2]]))
Int=np*M.determinant()
view(Int)

\newcommand{\Bold}[1]{\mathbf{#1}}1

---

Step 6: The program computes the integral (flux). Changes with problem - we must put in our intervals of integration.

integral(integral(Int,(v,0,-u+1)),(u,0,1))

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2}

---

---

So our flux is: $\iint_S x dydz+y \,dxdz+z \,dxdy = \frac{1}{2}=0.5$  

Notice that this last "makes sense" in our particular problem since the integrand is 1 and the area of the projected triangle is 1/2.

---

---

Scroll down to the bottom graph 2. It is a visualization of calculating the flux!

Graph 1: We must usually use the parameterization of S that we get to graph the surface. We can do this easily if the intervals for u and v are numbers.  

Here: $v \in [0,-u+1]$ and $u \in [0,1]$ so not numbers. I don't know if one can graph with intervals with variables, so I just "drew" the triangle a couple of different ways.

S=polygon([(1,0,0),(0,1,0),(0,0,1)])
show(S,figsize=3)

T=sum([line3d([(0,0,1),(c,1-c,0)], color=hue((c+8)/8), width=0.8) for c in [0..1,step=0.01]])
show(T,figsize=3)

---

Graph 2: Now we want a visual interpretation of what we have done. Remember that with the change of variables from the parametrization, we have:

• $\vec F= \lt u,v,1-u-v \gt$
• $\vec n=\lt 1,1,1 \gt$
• Integrand = 1

So the vectors that describe the normalized flux volume (and scaled by vector-scale=$vs$) are: 

• Starting point: $\lt u,v,1-u-v \gt$
• End point:  ${\lt u,v,1-u-v \gt }$ + $\lt \color{teal}{1} \cdot \color{purple}{1} \cdot vs,\color{teal}{1} \cdot \color{purple}{1}\cdot vs,\color{teal}{1} \cdot \color{purple}{1}\cdot vs \gt$= ${\lt u,v,1-u-v \gt }$ =${\lt u+1\cdot vs,v+1\cdot vs,1-u-v+1\cdot vs \gt }$

Remember that our surface with respect to u and v has $u \in [0,1]$ and $v \in [0,-u+1].

To do our programming interation, we substitute "c" for u and "d" for v and run d from [0,-c+1] and c from [0,1].  (Usually start tries with big steps and see how it goes and then reduce step size.)

The FLUX is the VOLUME between the surface S and the surface 'formed' by the endpoints of the arrows (multiplied by the orientation and divided by $\left\| {\vec n} \right\|$ )!

Here I have added a "blank" triangle so that the box is nice. Probably and easier way to do this...

vs=.5
vf=sum([sum([arrow3d((c,d,1-c-d),(c+1*vs,d+1*vs,1-c-d+1*vs), color=hue((c+8)/8), width=0.8) for d in [0..(-c+1),step=0.2]]) for c in [0..1,step=0.2]])
T_blank=polygon([(2,0,0),(0,2,0),(0,0,2)], color='white', opacity=0)
show(T+vf+T_blank,figsize=3)