Surface Integral of Flux Type - SOLVER

117 days ago by lfahlberg

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)
\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)
\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)
\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
\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)
\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
\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}
\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)