Integrating Functions

Indefinite integrals

Mathematica

Mathematica’s Integrate function can return the indefinite integral (antiderivative) of a function. Simply provide the function and the variable of integration.

In[1]:= Integrate[ x^3, x]

         4
        x
Out[1]= --
        4

Maxima

Maxima’s integrate() function can return the indefinite integral (antiderivative) of a function. Simply provide the function and the variable of integration.

(C1) integrate(x^2,x);

                                       3
                                      x
(D1)                                  --
                                      3

Axiom

Axiom’s function can return the indefinite integral (antiderivative) of a function. Simply provide the function and the variable of integration.

(1) -> integrate(x^2,x)
(1) ->
        1  3
   (1)  - x
        3
                                            Type: Polynomial Fraction Integer

Definite integrals

Mathematica

Mathematica’s Integrate function can return the definite integral of a function. The limits may be numbers or expressions. Simply provide the function together with the variable of integration and limits enclosed in a list.

In[1]:= Integrate[ x^2, {x,0,1}]

        1
Out[1]= -
        3

In[2]:= Integrate[ x^2,{x,y,y^2}]

          3    6
        -y    y
Out[2]= --- + --
         3    3

Maxima

Maxima’s integrate() function can return the definite integral of a function. The limits may be numbers or expressions. Simply provide the function, the variable of integration, and limits.

(C1) integrate( x^2, x, 0, 1);

                                       1
(D1)                                   -
                                       3
(C2) integrate( x^2, x, y, y^2);

                                     6    3
                                    y    y
(D2)                                -- - --
                                    3    3

Axiom

Axiom’s integrate() function can return the definite integral of a function. The limits may be numbers or expressions. Simply provide the function and the variable of integration with its range.

(1) -> integrate( x^2, x=0..1 )
(1) ->
        1
   (1)  -
        3
                    Type: Union(f1: OrderedCompletion Expression Integer,...)
(2) -> integrate( x^2, x=y..y^2 )
(2) ->
         6    3
        y  - y
   (2)  -------
           3
                    Type: Union(f1: OrderedCompletion Expression Integer,...)

Numerically computed integrals

Mathematica

Integrals may be computed numerically in Mathematica, by substituting NIntegrate[] for Integrate[] in definite integrals:

In[1]:= NIntegrate[E^(-x^2),{x,0,1}]

Out[1]= 0.746824

In[2]:= NIntegrate[  Sqrt[1+x^2+y^2], {y,0,1}, {x,0,1} ]

Out[2]= 1.28079

Maxima

Integrals may be computed numerically in Maxima using an algorithm called Romberg’s method. The syntax is precisely the same as for symbolically computed definite integrals, only substituting the function name romberg() for integrate():

(%i1) romberg( %e^(-x^2),x,0,1 );

(%o1)                           0.7468241699099
(%i2) romberg( romberg(sqrt(1+x^2+y^2),y,0,1), x,0,1 );

(%o2)                          1.280788828222111

There are other numerical integration methods available to Maxima. Another one is Newton-Coates integration. This can be important because some functions integrate better with different methods (it never hurts to try two methods and compare). The Newton-Coates method has to be loaded before used, via load(qq) and the function name is quanc8().

(%i1) load("qq");

(%o1)            /usr/share/maxima/5.9.1/share/numeric/qq.lisp
(%i2) quanc8( %e^(-x^2),x,0,1 );

(%o2)                          0.74682413281243
(%i3) quanc8( quanc8( sqrt(1+x^2+y^2),y,0,1), x,0,1 );

(%o3)                          1.280789275273411

Notice that it gives a result slightly different from romberg() for each integral.

Axiom

Axiom has several functions for computing integrals numerically, although I have only been able to find algorithms for computing single integrals. Unfortunately, all of the functions differ substantially from the usual syntax of integrate(), so they feel a bit foreign.

One function is romberg(), and the arguments in order are

the function (with its independent variable explicit)
the lower limit
the upper limit
a relative error bound
an absolute error bound
a minimum number of steps to compute
a maximum number of steps to compute

(1) -> romberg( x +-> %e^(-x^2), 0, 1, 1.e-6, 1.e-6, 6, 10 )
(1) ->
   (1)
   [value= 0.7468241328 124270254, error= 0.324 E -16, totalpts= 129,
    success= true]
   Type: Record(value: Float,error: Float,totalpts: Integer,success: Boolean)

The function has a success value of true if the answer is returned within either the relative or absolute error before the maximum number of steps is reached.

Octave

Octave doesn’t contain a native facility for doing numerical integration, but it is a simple matter to write a program that approximates an integral.

octave:1> function val = midpoint(f,a,b,n)
> # midpoint(f,a,b,n) approximate the integral of f over the
> # interval [a,b] using n rectangles.  f is the name of
> # the function to evaluate.
>
> h = (b-a)/n;
> x = [1:n]*h + a - h/2;
> val = sum(feval(f,x)) * h;
> endfunction
octave:2> function y = f(x); y = x.^2; endfunction
octave:3> midpoint("f", 0, 1, 50 )
ans = 0.33330

Double integrals in polar coordinates

Mathematica

I find it most convenient to do double integrals in polar coordinates if I get my steps in the correct order. Generally, defining x and y up front works best. Lets find the volume of the "ice cream" cone above the cone z = sqrt(x^2 + y^2) and below the sphere z = sqrt(1-x^2-y^2).

In[1]:= x = r*Cos[t]; y = r*Sin[t]

Out[1]= r Sin[t]

In[2]:= z1 = Simplify[ Sqrt[x^2+y^2]]

              2
Out[2]= Sqrt[r ]

In[3]:= z2 = Simplify[ Sqrt[1-x^2-y^2]]

                  2
Out[3]= Sqrt[1 - r ]

In[4]:= Solve[ z1==z2, r ]

Solve::ifun: Inverse functions are being used by Solve, so some solutions may
     not be found; use Reduce for complete solution information.

                    1                1
Out[4]= {{r -> -(-------)}, {r -> -------}}
                 Sqrt[2]          Sqrt[2]

Let’s extract just the 2nd solution, since that is the one we want for the upper limit for r, and integrate.

In[5]:= r /. Out[4][[2]]

           1
Out[5]= -------
        Sqrt[2]

In[6]:= Integrate[ (z2-z1)*r, {r,0,Out[5]} ]

        2 - Sqrt[2]
Out[6]= -----------
             6

Everything looks good, so we compute the outer integral, and we’re finished.

In[7]:= Integrate[ Out[6], {t,0,2*Pi} ]

        2 Pi   Sqrt[2] Pi
Out[7]= ---- - ----------
         3         3

Maxima

(%i1) x : r*cos(t); y : r*sin(t);
(%o1)                              r cos(t)
(%o2)                              r sin(t)
(%i3) z1 : trigsimp( sqrt(x^2+y^2));
(%o3)                                  r
(%i4) z2 : trigsimp( sqrt(1-x^2-y^2));
                                           2
(%o4)                            sqrt(1 - r )
(%i5) solve( z1=z2, r );
                                             2
(%o5)                         [r = sqrt(1 - r )]

At this step, Maxima is not being too helpful. We’re going to have to help it along or figure out what r should be on our own. If we square both sides, that is enough to nudge Maxima.

(%i6) solve( r^2 = 1-r^2, r );
                                   1            1
(%o6)                    [r = - -------, r = -------]
                                sqrt(2)      sqrt(2)

That’s better. Let’s extract just the 2nd solution, since that is the one we want for the upper limit for r, and integrate.

(%i7) r, %o6[2];
                                       1
(%o7)                               -------
                                    sqrt(2)
(%i8) integrate( (z2-z1)*r, r, 0, %o7 );
                                  1   sqrt(2)
(%o8)                             - - -------
                                  3      6

Everything looks good, so we compute the outer integral, and we’re finished.

(%i9) integrate( %o8, t, 0, 2*%pi );
                                 1   sqrt(2)
(%o9)                         2 (- - -------) %pi
                                 3      6

Curve integrals

Mathematica

Mathematica can do most of the tedious work involved in computing curve integrals and help prevent algebraic errors. It is convenient to use the Mag[] funtion listed on Vectors.html, and it is copied below.

In[1]:= Mag[a_] := Sqrt[ a . a ]

To calculate \int_C 2x \, ds where C is the part of the parabola y=x^2 that goes from (0,0) to (1,1), you might do:

In[2]:= x = t; y = t^2;

In[3]:= r = { x, y }

             2
Out[3]= {t, t }

In[4]:= dr = D[r,t]

Out[4]= {1, 2 t}

In[5]:= ds = Mag[dr]

                    2
Out[5]= Sqrt[1 + 4 t ]

In[6]:= Integrate[ 2*x*ds, {t,0,1}]

        -1 + 5 Sqrt[5]
Out[6]= --------------
              6

In[7]:= N[%]

Out[7]= 1.69672

Here’s a more elaborate example, for a helix in R³. Evaluate \int_C y \sin z \, ds where C is the helix < \cos t, \sin t, t> and 0 <= t <= 2*pi:

In[1]:= Mag[a_] := Sqrt[ a . a ]

In[2]:= x = Cos[t]; y = Sin[t]; z = t;

In[3]:= r = {x,y,z}

Out[3]= {Cos[t], Sin[t], t}

In[4]:= dr = D[r,t]

Out[4]= {-Sin[t], Cos[t], 1}

In[5]:= ds = Mag[dr]

                       2         2
Out[5]= Sqrt[1 + Cos[t]  + Sin[t] ]

You might note at this moment that ds can be simplified. You can ask Mathematica to perform simplifications with the Simplify[] function.

In[6]:= Simplify[%]

Out[6]= Sqrt[2]

Now redefine ds to be the simplified expression and carry on with the integral.

In[7]:= ds = %

Out[7]= Sqrt[2]

In[8]:= Integrate[ y*Sin[z]*ds, {t,0,2*Pi} ]

Out[8]= Sqrt[2] Pi

Maxima

Maxima can do most of the tedious work involved in computing curve integrals and help prevent algebraic errors. It is convenient to use the mag() function listed on Vectors.html, and it is copied below.

(%i1) mag(a) := sqrt( a . a );

(%o1)                        mag(a) := SQRT(a . a)

To calculate \int_C 2x \, ds where C is the part of the parabola y=x^2 that goes from (0,0) to (1,1), you might do:

(%i2) x : t; y : t^2;

(%o2)                                  t
(%i3)
                                       2
(%o3)                                 t
(%i4) r : [ x, y ];

                                         2
(%o4)                               [t, t ]
(%i5) dr : diff( r, t );

(%o5)                              [1, 2 t]
(%i6) ds : mag(dr);

                                        2
(%o6)                           SQRT(4 t  + 1)
(%i7) integrate( 2*x*ds, t, 0, 1 );

                                 5 SQRT(5)   1
(%o7)                         2 (--------- - --)
                                    12       12
(%i8) %,numer;

(%o8)                          1.696723314583158

Here’s a more elaborate example, for a helix in R³. Evaluate \int_C y \sin z \, ds where C is the helix < \cos t, \sin t, t> and 0 <= t <= 2*pi:

(%i1) mag(a) := sqrt( a . a );

(%o1)                        mag(a) := SQRT(a . a)
(%i2) x : cos(t); y : sin(t); z : t;

(%o2)                               COS(t)
(%i3)
(%o3)                               SIN(t)
(%i4)
(%o4)                                  t
(%i5) r : [ x, y, z ];

(%o5)                         [COS(t), SIN(t), t]
(%i6) dr : diff(r,t);

(%o6)                        [- SIN(t), COS(t), 1]
(%i7) ds : mag(dr);

                                  2         2
(%o7)                     SQRT(SIN (t) + COS (t) + 1)

You might note at this moment that ds can be simplified. You can ask Maxima to perform trigonometric simplifications with the trigsimp() function.

(%i8) trigsimp(ds);

(%o8)                               SQRT(2)

Now redefine ds to be the simplified expression and carry on with the integral.

(%i9) ds : %;

(%o9)                               SQRT(2)
(%i10) integrate( y*sin(z)*ds, t, 0, 2*%pi );

(%o10)                            SQRT(2) %PI

Axiom

Axiom can do most of the tedious work involved in computing curve integrals and help prevent algebraic errors. However, one ability that Axiom lacks is the ability to conveniently take the derivative of a vector (and have the derivative apply to each component). So a new DV() function is introduced for that:

(1) -> macro DV(vec,var) == map( f +-> D(f,var), vec )
                                                                   Type: Void

To calculate \int_C 2x \, ds where C is the part of the parabola y=x^2 that goes from (0,0) to (1,1), you might do:

(2) -> x := t; y := t^2
(2) ->
         2
   (2)  t
                                                     Type: Polynomial Integer
(3) -> r := vector( [x,y] )
(3) ->
            2
   (3)  [t,t ]
                                              Type: Vector Polynomial Integer
(4) -> dr := DV(r,t)
(4) ->
   (4)  [1,2t]
                                              Type: Vector Polynomial Integer
(5) -> ds := length(dr)
(5) ->
         +-------+
         |  2
   (5)  \|4t  + 1
                                                     Type: Expression Integer
(6) -> integrate( 2*x*ds, t=0..1 )
(6) ->
   (6)  potentialPole
                                         Type: Union(pole: potentialPole,...)

Axiom is afraid, because of the square root, that there might be a pole (a place where the antiderivative goes to infinity). We need to reassure Axiom that everything is all right to integrate.

(7) -> integrate( 2*x*ds, t=0..1, "noPole" )
(7) ->
              +-+
        - 207\|5  + 463
   (7)  ---------------
             +-+
         102\|5  - 228
                    Type: Union(f1: OrderedCompletion Expression Integer,...)
(8) -> % :: Expression Float
(8) ->
   (8)  1.6967233145 831580375
                                                       Type: Expression Float

Here’s a more elaborate example, for a helix in R³. Evaluate \int_C y \sin z \, ds where C is the helix < \cos t, \sin t, t > and 0 <= t <= 2*pi:

(1) -> macro DV(vec,var) == map( f +-> D(f,var), vec )
                                                                   Type: Void
(2) -> x := cos(t); y := sin(t); z := t;
(2) ->
                                                             Type: Variable t
(3) -> r := vector( [x,y,z] )
(3) ->
   (3)  [cos(t),sin(t),t]
                                              Type: Vector Expression Integer
(4) -> dr := DV(r,t)
(4) ->
   (4)  [- sin(t),cos(t),1]
                                              Type: Vector Expression Integer
(5) -> ds := length(dr)
(5) ->
         +---------------------+
         |      2         2
   (5)  \|sin(t)  + cos(t)  + 1
                                                     Type: Expression Integer

You might note at this moment that ds can be simplified. You can ask Axiom to perform simplifications with the simplify() function.

(6) -> simplify(ds)
(6) ->
         +-+
   (6)  \|2
                                                     Type: Expression Integer

Now redefine ds to be the simplified expression and carry on with the integral.

(7) -> ds := %
(7) ->
         +-+
   (7)  \|2
                                                     Type: Expression Integer
(8) -> integrate( y*sin(z)*ds, t=0..2*%pi, "noPole" )
(8) ->
            +-+
   (8)  %pi\|2
                    Type: Union(f1: OrderedCompletion Expression Integer,...)

Work done by a vector field

Mathematica

The work done by a vector field around a curve is a special form of curve integral, and Mathematica can help calculate it. For example, let us evaluate the line integral \int_C F \cdot dr, where F(x,y) = < x^2*y^3, -y*\sqrt(x) > and C is parametrized by r(t) = < t^2, t^3 > on the interval [0,1].

In[1]:= x = t^2; y = t^3;

In[2]:= r = {x,y}

          2   3
Out[2]= {t , t }

In[3]:= dr = D[r,t]

                 2
Out[3]= {2 t, 3 t }

In[4]:= F = { x^2*y^3, -y*Sqrt[x] }

          13     3       2
Out[4]= {t  , -(t  Sqrt[t ])}

In[5]:= Integrate[ F . dr, {t,0,1} ]

          31
Out[5]= -(---)
          105

Maxima

The work done by a vector field around a curve is a special form of curve integral, and Maxima can help calculate it. For example, let us evaluate the line integral \int_C F \cdot dr, where F(x,y) = < x^2*y^3, -y*\sqrt(x) > and C is parametrized by r(t) = < t^2, t^3 > on the interval [0,1].

(%i1) x : t^2; y : t^3;

                                       2
(%o1)                                 t
(%i2)
                                       3
(%o2)                                 t
(%i3) r : [x,y];

                                     2   3
(%o3)                              [t , t ]
(%i4) dr : diff(r,t);

                                           2
(%o4)                             [2 t, 3 t ]
(%i5) F : [ x^2*y^3, -y*sqrt(x) ];

                                13     3
(%o5)                         [t  , - t  ABS(t)]
(%i6) integrate( F . dr, t, 0, 1 );

                                       31
(%o6)                                - ---
                                       105

Axiom

(1) -> macro DV(vec,var) == map( f +-> D(f,var), vec )
                                                                   Type: Void
(2) -> x := t^2; y := t^3;
(2) ->
                                                     Type: Polynomial Integer
(3) -> r := vector [x,y]
(3) ->
          2  3
   (3)  [t ,t ]
                                              Type: Vector Polynomial Integer
(4) -> dr := DV(r,t)
(4) ->
              2
   (4)  [2t,3t ]
                                              Type: Vector Polynomial Integer
(5) -> F := vector( [ x^2*y^3, -y*sqrt(x) ] )
(5) ->
                  +--+
          13    3 | 2
   (5)  [t  ,- t \|t  ]
                                              Type: Vector Expression Integer
(6) -> integrate( dot( F, dr ), t=0..1, "noPole" )
(6) ->
           31
   (6)  - ---
          105
                    Type: Union(f1: OrderedCompletion Expression Integer,...)

Surface integrals

Mathematica

Mathematica can do most of the tedious work involved in computing surface integrals and help prevent algebraic errors.

Calculate \int_S x^2 \, ds where S is the sphere of radius 1, centered at the origin. This surface is most naturally parametrized in terms of spherical coordinates.

In[1]:= x = Cos[t]*Sin[p]; y = Sin[t]*Sin[p]; z = Cos[p]

Out[1]= Cos[p]

In[2]:= r = { x,y,z }

Out[2]= {Cos[t] Sin[p], Sin[p] Sin[t], Cos[p]}

In[3]:= rt = D[r,t]

Out[3]= {-(Sin[p] Sin[t]), Cos[t] Sin[p], 0}

In[4]:= rp = D[r,p]

Out[4]= {Cos[p] Cos[t], Cos[p] Sin[t], -Sin[p]}

In[5]:= dS = Norm[ Cross[ rt, rp ]]

                              2 2             2        2
Out[5]= Sqrt[Abs[Cos[t] Sin[p] ]  + Abs[Sin[p]  Sin[t]]  +

                         2                               2 2
>     Abs[-(Cos[p] Cos[t]  Sin[p]) - Cos[p] Sin[p] Sin[t] ] ]

In[6]:= Integrate[ x^2*dS, {t,0,2*Pi} ]

                                                             2        2
        Pi Sqrt[-(Cos[2 Re[p]] Cosh[2 Im[p]]) + Cosh[2 Im[p]] ] Sin[p]
Out[6]= ---------------------------------------------------------------
                                    Sqrt[2]

In[7]:= Integrate[ %, {p,0,Pi}]

        4 Pi
Out[7]= ----
         3

Maxima

Maxima can do most of the tedious work involved in computing surface integrals and help prevent algebraic errors. It is convenient to use the cross() and mag() functions listed on Vectors.html. They are copied below.

(%i1) cross(a,b) := [ a[2]*b[3]-a[3]*b[2], a[3]*b[1]-a[1]*b[3], a[1]*b[2]-a[2]*b[1] ];

(%o1)    cross(a, b) := [a  b  - a  b , a  b  - a  b , a  b  - a  b ]
                          2  3    3  2   3  1    1  3   1  2    2  1
(%i2) mag(a) := sqrt( a . a );

(%o2)                        mag(a) := SQRT(a . a)

Calculate \int_S x^2 \, ds where S is the sphere of radius 1, centered at the origin. This surface is most naturally parametrized in terms of spherical coordinates.

(%i3) x : cos(t)*sin(p); y : sin(t)*sin(p); z : cos(p);

(%o3)                            SIN(p) COS(t)
(%i4)
(%o4)                            SIN(p) SIN(t)
(%i5)
(%o5)                               COS(p)
(%i6) r : [ x, y, z ];

(%o6)               [SIN(p) COS(t), SIN(p) SIN(t), COS(p)]
(%i7) rt : diff(r,t);

(%o7)                 [- SIN(p) SIN(t), SIN(p) COS(t), 0]
(%i8) rp : diff(r,p);

(%o8)              [COS(p) COS(t), COS(p) SIN(t), - SIN(p)]
(%i9) dS : mag(cross(rt,rp));

                               2                       2    2      4       2
(%o9) SQRT((- COS(p) SIN(p) SIN (t) - COS(p) SIN(p) COS (t))  + SIN (p) SIN (t)

                                                                  4       2
                                                             + SIN (p) COS (t))

It turns out that dS is a bit of a mess, but simplifies. You can get Maxima to do some basic trig simplifying by applying the trigsimp() function.

(%i10) trigsimp(%);

(%o10)                              SIN(p)

Redefine dS as this simplified expression and we’re ready to integrate with respect to each variable in turn.

(%i11) dS : %;
(%o11)                              SIN(p)
(%i12) integrate( x^2*dS, t, 0, 2*%pi );

                                         3
(%o12)                            %PI SIN (p)
(%i13) integrate( %, p, 0, %pi );

                                     4 %PI
(%o13)                               -----
                                       3

Axiom

Axiom can do most of the tedious work involved in computing surface integrals and help prevent algebraic errors. It uses the DV() macro listed on Vectors.html, copied below.

(1) -> macro DV(vec,var) == map( f +-> D(f,var), vec )
                                                                   Type: Void

Calculate \int_S x^2 \, ds where S is the sphere of radius 1, centered at the origin. This surface is most naturally parametrized in terms of spherical coordinates.

(2) -> x := cos(t)*sin(p); y := sin(t)*sin(p); z := cos(p)
(2) ->
   (2)  cos(p)
                                                     Type: Expression Integer
(3) -> r := vector [x,y,z]
(3) ->
   (3)  [cos(t)sin(p),sin(p)sin(t),cos(p)]
                                              Type: Vector Expression Integer
(4) -> rt := DV(r,t)
(4) ->
   (4)  [- sin(p)sin(t),cos(t)sin(p),0]
                                              Type: Vector Expression Integer
(5) -> rp := DV(r,p)
(5) ->
   (5)  [cos(p)cos(t),cos(p)sin(t),- sin(p)]
                                              Type: Vector Expression Integer
(6) -> dS := length( cross(rt,rp))
(6) ->
   (6)
   ROOT
              2      2      4          4          2      2      2       2
        cos(p) sin(p) sin(t)  + (sin(p)  + 2cos(p) cos(t) sin(p) )sin(t)
      +
              2      4         2      4      2
        cos(t) sin(p)  + cos(p) cos(t) sin(p)
                                                     Type: Expression Integer

It turns out that dS is a bit of a mess, but simplifies. You can get Axiom to do some simplifying by applying the simplify() function. Then we will be ready to integrate with respect to each variable in turn.

(7) -> dS := simplify( % )
(7) ->
         +-------------+
         |        2
   (7)  \|- cos(p)  + 1
                                                     Type: Expression Integer
(8) -> integrate( x^2*dS, t=0..2*%pi )
(8) ->
   (8)  potentialPole
                                         Type: Union(pole: potentialPole,...)

Because of the square root, Axiom thinks there might be a pole (a place where the antiderivative has an asymptote). Recompute with the "noPole" argument.

(9) -> integrate( x^2*dS, t=0..2*%pi, "noPole" )
(9) ->
                    +-------------+
                  2 |        2
   (9)  %pi sin(p) \|- cos(p)  + 1
                    Type: Union(f1: OrderedCompletion Expression Integer,...)
(10) -> integrate( %, p=0..%pi, "noPole" )
(10) ->
        4%pi
   (10)  ----
          3
                    Type: Union(f1: OrderedCompletion Expression Integer,...)

Flux integral

Mathematica

Really a special kind of surface integral, a flux integral can be calculated with Mathematica. Let F(x,y,z) = < 0, 1-x^2, 0 > be the velocity at which a fluid flows. What is the rate it which the fluid flows through the unit disk x^2+z^2=1? It makes sense to parametrize with polar coordinates (u=radius, v=angle).

In[1]:= x = u*Cos[v]; y = 0; z = u*Sin[v]

Out[1]= u Sin[v]

In[2]:= r = {x,y,z}

Out[2]= {u Cos[v], 0, u Sin[v]}

In[3]:= ru = D[r,u]

Out[3]= {Cos[v], 0, Sin[v]}

In[4]:= rv = D[r,v]

Out[4]= {-(u Sin[v]), 0, u Cos[v]}

In[5]:= dS = Cross[ru,rv]

                      2            2
Out[5]= {0, -(u Cos[v] ) - u Sin[v] , 0}

Since this could use a bit of simplification, we’ll take care of that before going on.

In[6]:= Simplify[%]

Out[6]= {0, -u, 0}

In[7]:= dS = %

Out[7]= {0, -u, 0}

Now to finish,

In[8]:= F = { 0, 1-x^2, 0 }

                 2       2
Out[8]= {0, 1 - u  Cos[v] , 0}

In[9]:= Integrate[ F . dS, {u,0,1} ]

                     2
          1    Cos[v]
Out[9]= -(-) + -------
          2       4

In[10]:= Integrate[ %, {v,0,2*Pi} ]

         -3 Pi
Out[10]= -----
           4

In this case, the negative is not significant, merely an indication that we oriented our circle against the flow of the fluid.

Maxima

Really a special kind of surface integral, a flux integral can be calculated with Maxima. We’ll want the cross product function from Vectors.html to evaluate this:

(%i1) cross(a,b) := [ a[2]*b[3]-a[3]*b[2], a[3]*b[1]-a[1]*b[3], a[1]*b[2]-a[2]*b[1] ];

(%o1)    cross(a, b) := [a  b  - a  b , a  b  - a  b , a  b  - a  b ]
                          2  3    3  2   3  1    1  3   1  2    2  1

Let F(x,y,z) = < 0, 1-x^2, 0 > be the velocity at which a fluid flows. What is the rate it which the fluid flows through the unit disk x^2+z^2=1? It makes sense to parametrize with polar coordinates (u=radius, v=angle).

(%i2) x : u*cos(v); y:0; z:u*sin(v);

(%o2)                              u COS(v)
(%i3)
(%o3)                                  0
(%i4)
(%o4)                              u SIN(v)
(%i5) r : [x,y,z];

(%o5)                       [u COS(v), 0, u SIN(v)]
(%i6) ru : diff(r,u);

(%o6)                         [COS(v), 0, SIN(v)]
(%i7) rv : diff(r,v);

(%o7)                      [- u SIN(v), 0, u COS(v)]
(%i8) dS : cross(ru,rv);

                                   2           2
(%o8)                   [0, - u SIN (v) - u COS (v), 0]

Since this could use a bit of simplification, we’ll take care of that before going on.

(%i9) trigsimp(%);

(%o9)                             [0, - u, 0]
(%i10) dS : %;

(%o10)                            [0, - u, 0]

Now to finish,

(%i11) F : [ 0, 1-x^2, 0 ];

                                     2    2
(%o11)                      [0, 1 - u  COS (v), 0]
(%i12) integrate( F . dS, u, 0, 1 );

                         4           2
                      COS (v) - 2 COS (v) + 1       1
(%o12)                ----------------------- - ---------
                                  2                  2
                             4 COS (v)          4 COS (v)
(%i13) integrate( %, v, 0, 2*%pi );

                                      3 %PI
(%o13)                              - -----
                                        4

In this case, the negative is not significant, merely an indication that we oriented our circle against the flow of the fluid.

Axiom

Really a special kind of surface integral, a flux integral can be calculated with Axiom. We’ll want the DV() macro from Vectors.html to evaluate this:

(1) -> macro DV(vec,var) == map( f +-> D(f,var), vec )
                                                                   Type: Void

(2) -> x := u*cos(v); y := 0; z := u*sin(v)
(2) ->
   (2)  u sin(v)
                                                     Type: Expression Integer
(3) -> r := vector [x,y,z]
(3) ->
   (3)  [u cos(v),0,u sin(v)]
                                              Type: Vector Expression Integer
(4) -> ru := DV(r,u)
(4) ->
   (4)  [cos(v),0,sin(v)]
                                              Type: Vector Expression Integer
(5) -> rv := DV(r,v)
(5) ->
   (5)  [- u sin(v),0,u cos(v)]
                                              Type: Vector Expression Integer
(6) -> dS := cross(ru,rv)
(6) ->
                     2           2
   (6)  [0,- u sin(v)  - u cos(v) ,0]
                                              Type: Vector Expression Integer

This could use a bit of simplification, but Axiom’s simplify() function doesn’t work on vectors. One way to proceed is to use the map() function to apply simplify() to each of the expressions in the vector.

(7) -> map( f +-> simplify(f), % )
(7) ->
   (7)  [0,- u,0]
                                              Type: Vector Expression Integer
(8) -> dS := %
(8) ->
   (8)  [0,- u,0]
                                              Type: Vector Expression Integer

Now to finish,

(9) ->  F := vector [ 0, 1-x^2, 0 ]
(9) ->
              2      2
   (9)  [0,- u cos(v)  + 1,0]
                                              Type: Vector Expression Integer
(10) -> integrate( dot(F,dS), u=0..1 )
(10) ->
               2
         cos(v)  - 2
   (10)  -----------
              4
                    Type: Union(f1: OrderedCompletion Expression Integer,...)
(11) -> integrate( %, v=0..2*%pi )
(11) ->
   (11)  potentialPole
                                         Type: Union(pole: potentialPole,...)

Assure Axiom that there is no pole.

(12) -> integrate( %%(10), v=0..2*%pi, "noPole" )
(12) ->
           3%pi
   (12)  - ----
             4
                    Type: Union(f1: OrderedCompletion Expression Integer,...)

In this case, the negative is not significant, merely an indication that we oriented our circle against the flow of the fluid.

Stokes' integral

Maxima

The double integral that occurs in Stokes' theorem can be calculated on the computer almost exactly like other Flux integrals, but there is one small catch that is noted below.

As our example, let ‘F(x,y,z) = < -y, x, xy >` be a vector field. We want to computer the work around the unit circle x^2 + y^2 = 1. To make use of Stokes’ theorem, we’ll compute a double integral over the top half of the unit sphere.

Like with other surface integrals, we’ll want the cross product function from Vectors.html to evaluate this:

(%i1) cross(a,b) := [ a[2]*b[3]-a[3]*b[2], a[3]*b[1]-a[1]*b[3], a[1]*b[2]-a[2]*b[1] ];

(%o1)    cross(a, b) := [a  b  - a  b , a  b  - a  b , a  b  - a  b ]
                          2  3    3  2   3  1    1  3   1  2    2  1

With surface integrals, it is usually simplest to start with a parameterization. In this case, the reasonable choice is spherical coordinates, phi and theta. We’ll use p and t respectively.

(%i2) x : cos(t)*sin(p); y : sin(t)*sin(p); z : cos(p);

(%o2)                            SIN(p) COS(t)
(%i3)
(%o3)                            SIN(p) SIN(t)
(%i4)
(%o4)                               COS(p)
(%i5) r : [x,y,z];

(%o5)               [SIN(p) COS(t), SIN(p) SIN(t), COS(p)]
(%i6) rt : diff(r,t); rp : diff(r,p); dS : cross(rt,rp);

(%o6)                 [- SIN(p) SIN(t), SIN(p) COS(t), 0]
(%i7)
(%o7)              [COS(p) COS(t), COS(p) SIN(t), - SIN(p)]
(%i8)
            2                 2
(%o8) [- SIN (p) COS(t), - SIN (p) SIN(t),

                                                  2                       2
                               - COS(p) SIN(p) SIN (t) - COS(p) SIN(p) COS (t)]
(%i9) dS : trigsimp(dS);

                   2                 2
(%o9)        [- SIN (p) COS(t), - SIN (p) SIN(t), - COS(p) SIN(p)]

Now, we need to compute the curl of the vector field. We could use the curl function from Vectors.html.

(%i10) curl(a) := [
  diff(a[3],'y)-diff(a[2],'z),
  diff(a[1],'z)-diff(a[3],'x),
  diff(a[2],'x)-diff(a[1],'y) ];

(%o10) curl(a) := [DIFF(a , 'y) - DIFF(a , 'z), DIFF(a , 'z) - DIFF(a , 'x),
                         3              2             1              3

                                                   DIFF(a , 'x) - DIFF(a , 'y)]
                                                         2              1

But there is a subtlety. If we make a naive definition of F in terms of x, y, z, Maxima will use the definitions above and translate directly everything to t and p and when we calculate the curl we will get zero (nothing is in terms of x or y or z so each derivative will be zero).

The solution to this dilemma is to let Maxima know we don’t want it to evaluate fully. The way this is done is to place a singe quote mark in front of each of the variables.

(%i11) F : [ -'y, 'x, 'x*'y ];

(%o11)                           [- y, x, x y]
(%i12) cF : curl(F);

(%o12)                            [x, - y, 2]

At this point, we can use cF directly in our integral, or first explicity ask Maxima to evaluate it.

(%i13) cF : cF,eval;

(%o13)                [SIN(p) COS(t), - SIN(p) SIN(t), 2]
(%i14) integrate( cF . dS, t, 0, 2*%pi );

(%o14)                       - 4 %PI COS(p) SIN(p)
(%i15) integrate( %, p, 0, %pi/2 );

(%o15)                              - 2 %PI

If you understand this well, you’ll note that our answer was actually the negative of what we should have gotten. The reason for this can be seen in line (%o9). We see that the normal vectors point to the origin, but the properly oriented surface would have normal vectors that point out (i.e. they sholdn’t have those minus signs). This is just a matter of managing our parameterization poorly. If we had computed instead cross(rp,rt) in line (%i6) the sign would have been correct.