Application Center - Maplesoft

App Preview:

The VectorCalculus Package

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

The VectorCalculus Package in Maple 8

Copyright 2002 Waterloo Maple Inc.

Maple 8 provides a new package called VectorCalculus  for computing with vectors, vector fields, multivariate functions and parametric curves.  Computations include vector arithmetic using basis vectors, "del"-operations, multiple integrals over regions and solids, line and surface integrals, differential-geometric properties of curves, and many others.  

You can easily perform these computations in any coordinate system and convert results between coordinate systems. The package is fully compatible with the Maple LinearAlgebra  package.  You can also extend the VectorCalculus  package by defining your own coordinate systems.

>    restart:
with(VectorCalculus);

Warning, the assigned names <,> and <|> now have a global binding

Warning, these protected names have been redefined and unprotected: *, +, ., Vector, diff, int, limit, series

[`&x`, `*`, `+`, `.`, `<,>`, `<|>`, AddCoordinates, ArcLength, BasisFormat, Binormal, CrossProduct, Curl, Curvature, Del, DirectionalDiff, Divergence, DotProduct, Flux, GetCoordinateParameters, GetCoor...
[`&x`, `*`, `+`, `.`, `<,>`, `<|>`, AddCoordinates, ArcLength, BasisFormat, Binormal, CrossProduct, Curl, Curvature, Del, DirectionalDiff, Divergence, DotProduct, Flux, GetCoordinateParameters, GetCoor...
[`&x`, `*`, `+`, `.`, `<,>`, `<|>`, AddCoordinates, ArcLength, BasisFormat, Binormal, CrossProduct, Curl, Curvature, Del, DirectionalDiff, Divergence, DotProduct, Flux, GetCoordinateParameters, GetCoor...
[`&x`, `*`, `+`, `.`, `<,>`, `<|>`, AddCoordinates, ArcLength, BasisFormat, Binormal, CrossProduct, Curl, Curvature, Del, DirectionalDiff, Divergence, DotProduct, Flux, GetCoordinateParameters, GetCoor...

>   

Creating Vectors and Vector Fields

It's usually convenient to specify a default coordinate system using the SetCoordinates  function.  (For a listing of the 37 coordinate systems predefined in Maple, see the help page ?Coordinates).

>    SetCoordinates( cartesian[x, y, z] );

cartesian[x,y,z]

You can define vectors in the same way you do with Maple's LinearAlgebra  package, using either the Vector  constructor, or <>  notation.  
Note the output appears in terms of the cartesian basis vectors
e[x] , e[y]  and e[z] .  Vectors defined this way are recognized by the LinearAlgebra  pakcage.

>    v := Vector( [a, b, c] );

v := Vector(%id = 19703984)

The MapToBasis  function converts vectors and vector fields to other coordinate systems.

>    MapToBasis( v, cylindrical[r, theta, z] );

Vector(%id = 20231752)

A vector field is a data type in the VectorCalculus  package. The bars over the basis vectors distinguish vector fields from constant vectors.

>    F := VectorField( <y, -x, 1> );

F := Vector(%id = 20260688)

>    attributes( F );

vectorfield, coords = cartesian[x,y,z]

Evaluate the vector field F  at a particular vector, say `<,>`(4,b, 1) , using the evalVF  function

>    evalVF( F, <4, b, 1> );

Vector(%id = 20225056)

Change the default coordinate system to polar.

>    SetCoordinates( polar[r, theta] );

polar[r,theta]

Define a vector field in the polar coordinate system.

>    F := VectorField( <0, 1> );

F := Vector(%id = 20225136)

Below we see that the value of the basis vector e[theta]  is not constant over the plane.  This is a major difference from cartesian coordinates, in which the basis vectors e[x]  and e[y]  are constants over the plane (i.e. `<,>`(1,0)  and `<,>`(0,1) .)  You can use Maple to clarify this common point of confusion for students.

>    MapToBasis( F, cartesian[x, y]);

Vector(%id = 20475980)

>   

Operations on Vectors and Vector Fields

Basic Arithmetic

The package is not limited to three dimensions.  The cartesian  option without variables specified allows you to compute with vectors of any dimension.

>    SetCoordinates( cartesian );

cartesian

Vector addition and scalar multiplication over R^4

>    v1 := a * <1, 2, 5, -1> + b * <3, 11, 0, 8>;

v1 := Vector(%id = 20590944)

>    v2 := 2*v1 - v1/3;

v2 := Vector(%id = 20910808)

>   

Dot Product, Divergence and Del

The .  operator and the DotProduct  function are equivalent ways to compute vector dot products.

>    <a, b, c, d> . <e, f, g, h>;

a*e+b*f+c*g+d*h

>    DotProduct( <a, b, c, d>, <e, f, g, h> );

a*e+b*f+c*g+d*h

You can compute the divergence of a vector field using either the Divergence  function, or the dot product with the Del  or Nabla  operators.

>    F := VectorField( <x^2, y^2, z^2>, cartesian[x,y,z] );

F := Vector(%id = 20591344)

>    Divergence( F );

2*x+2*y+2*z

>    Del . F;

2*x+2*y+2*z

>    Nabla . F;

2*x+2*y+2*z

>   

Cross Product, Curl and Del

To compute the cross product of two vectors in R^3 , use the CrossProduct  function or the &x  operator.

>    CrossProduct( <a, b, c>, <d, e, f>  );

Vector(%id = 21099068)

>    <a, b, c> &x <d,e,f> ;

Vector(%id = 21148040)

Compute the curl of a vector field F  using either the Curl  function or the expression Del &x F.

Note that Maple distinguishes the vector  0  from the number 0 by writing a basis-vector next to the 0.

>    Curl( VectorField( <x^2, y^2, z^2>, cartesian[x,y,z] ) );

Vector(%id = 21077824)

Here's a nonconservative (rotational) vector field.

>    F := VectorField( <y, -x, z>, cartesian[x,y,z] );

F := Vector(%id = 21077904)

We know that it's nonconservative because its curl is nonzero.

>    Del &x F;

Vector(%id = 21077944)

Routines from the plots  package can accept outputs from the VectorCalculus  package.

>    plots[fieldplot3d]( F, x=-1..1, y=-1..1, z=-1/2..1/2, scaling=constrained, shading=zhue );

[Maple Plot]

Multivariate Differential Calculus

Earlier versions of Maple could already perform multivariate differential calculus, but the VectorCalculus  package makes it especially simple and intuitive.  Working with gradients, laplacians and hessians of real-valued functions is straightforward, as the following worked examples show.

Laplacians, Gradients, and Hessians, Oh My!

>    SetCoordinates( cartesian ):

Example 1: Find the unit vector normal to the surface z = x^2*y^2+y+1  at the point p = `<,>`(0,0,1)

We'll solve this problem as students are taught.  Introduce the 3-variable function f(x,y,z) = -z+x^2*y^2+y+1  and treat the surface as the level surface f(x,y,z) = 0

>    f := -z + x^2*y^2 + y + 1;

f := -z+x^2*y^2+y+1

>    p := <0,0,1>;

p := Vector(%id = 21944780)

The normal to the surface at p  is just (grad f)(p) , since the gradient of f  is normal to f' s level surfaces.  We compute the gradient with the function Gradient ...

>    Gf := Gradient( f, cartesian[x,y,z]);

Gf := Vector(%id = 21944940)

...and evaluate the gradient at p  with the evalVF  function.

>    Gf0 := evalVF( Gf, p );

Gf0 := Vector(%id = 21945100)

Finally, we normalize the gradient vector with the Normalize  function from the LinearAlgebra  package.

>    Gf0_unit := LinearAlgebra :- Normalize(Gf0, 2);

Gf0_unit := Vector(%id = 2443196)

A plot of the result.

>    a := plots[arrow]( p, Gf0_unit, color=blue, thickness=3):
S := plot3d( f+z, x=-1..1, y=-1..1, style=wireframe, axes=boxed):
plots[display](S, a, scaling=constrained, orientation=[-20,65]);

[Maple Plot]

>   

Example 2: Find the point on the paraboloid z = (x^2+y^2)/5  closest to the point p = `<,>`(1,1,0)

>    SetCoordinates(cartesian[x,y,z]):


We'll solve this as a constrained minimization problem using the Lagrange Multiplier Theorem.  
Interpret the paraboloid as a level surface of the function
g(x,y,z) = (x^2+y^2)/5-z .

>    g := (x^2+y^2)/5-z;

g := 1/5*x^2+1/5*y^2-z

>    p := <1,1,0>;

p := Vector(%id = 2626880)

Let the objective function f  be the squared distance from point p  to a point `<,>`(x,y, z) .  We wish to minimize f .

>    f := (x-p[1])^2 + (y-p[2])^2 + (z-p[3])^2;

f := (x-1)^2+(y-1)^2+z^2

Form the Lagrangian for the constrained minimization problem, f(x,y,z)-lambda*g(x,y,z)

>    L := f - lambda * g;

L := (x-1)^2+(y-1)^2+z^2-lambda*(1/5*x^2+1/5*y^2-z)

Compute the gradient of the Lagrangian.

>    GL := Gradient( L );

GL := Vector(%id = 2627080)

Solve the four equations g = 0, GL = 0   for the four unknowns x, y, z   and   lambda

>    sol := fsolve( {g = 0, seq( GL[i] = 0, i=1..3 )},
               {x,y,z,lambda} ) ;

sol := {lambda = -.6307977539, x = .8879736440, y = .8879736440, z = .3153988770}

Verify that the solution is a minimum using the generalized 2nd-derivative test, i.e., det(Hessian(f-lambda*g)) < 0

>    HL := Hessian( L, [x,y,z,lambda]);

HL := Matrix(%id = 19876476)

>    LinearAlgebra :- Determinant(HL);

-4+8/5*lambda-16/25*y^2-4/25*lambda^2+16/125*lambda*y^2-16/25*x^2+16/125*x^2*lambda

The determinant is negative at the solution found, so the solution is a minimum.

>    eval( %, sol );

-6.209547599

A plot of the result.

>    P := plot3d( g+z, x=0..2, y=0..2, style=wireframe, axes=boxed ):
L := plottools[line]( [1,1,0], eval([x,y,z], sol), thickness=3, color=red):
plots[display]( P, L, scaling=constrained, orientation=[30,85]);

[Maple Plot]

>   

Example 3: Solve the time-dependent heat equation for a cylindrical heat source with external heat generation

We find a general solution to the time-dependent heat equation
[Maple OLE 2.0 Object]

in cylindrical coordinates, where
u(r,theta,z,t)  is the temperature, q   is the thermal source strength, and k  is the thermal diffusivity.

The Laplacian  function is convenient for writing down the PDEs of classical physics, especially in non-cartesian coordinate systems.

>    Lap := Laplacian( u(r, theta, z, t), cylindrical[r,theta,z]);

Lap := 1/r*(diff(u(r,theta,z,t),r)+r*diff(u(r,theta,z,t),`$`(r,2))+1/r*diff(u(r,theta,z,t),`$`(theta,2))+r*diff(u(r,theta,z,t),`$`(z,2)))

>    HeatEqn := Lap + q/k = diff(u(r,theta,z, t), t);

HeatEqn := 1/r*(diff(u(r,theta,z,t),r)+r*diff(u(r,theta,z,t),`$`(r,2))+1/r*diff(u(r,theta,z,t),`$`(theta,2))+r*diff(u(r,theta,z,t),`$`(z,2)))+q/k = diff(u(r,theta,z,t),t)

Solve the partial differential equation.

>    pdsolve( HeatEqn, build );

u(r,theta,z,t) = exp(_c[3]^(1/2)*z)*exp(_c[2]^(1/2)*theta)*_C3*exp(_c[4]*t)*_C7*_C5*_C1*BesselJ((-_c[2])^(1/2),(_c[3]-_c[4])^(1/2)*r)+exp(_c[3]^(1/2)*z)*exp(_c[2]^(1/2)*theta)*_C3*exp(_c[4]*t)*_C7*_C5*...
u(r,theta,z,t) = exp(_c[3]^(1/2)*z)*exp(_c[2]^(1/2)*theta)*_C3*exp(_c[4]*t)*_C7*_C5*_C1*BesselJ((-_c[2])^(1/2),(_c[3]-_c[4])^(1/2)*r)+exp(_c[3]^(1/2)*z)*exp(_c[2]^(1/2)*theta)*_C3*exp(_c[4]*t)*_C7*_C5*...
u(r,theta,z,t) = exp(_c[3]^(1/2)*z)*exp(_c[2]^(1/2)*theta)*_C3*exp(_c[4]*t)*_C7*_C5*_C1*BesselJ((-_c[2])^(1/2),(_c[3]-_c[4])^(1/2)*r)+exp(_c[3]^(1/2)*z)*exp(_c[2]^(1/2)*theta)*_C3*exp(_c[4]*t)*_C7*_C5*...
u(r,theta,z,t) = exp(_c[3]^(1/2)*z)*exp(_c[2]^(1/2)*theta)*_C3*exp(_c[4]*t)*_C7*_C5*_C1*BesselJ((-_c[2])^(1/2),(_c[3]-_c[4])^(1/2)*r)+exp(_c[3]^(1/2)*z)*exp(_c[2]^(1/2)*theta)*_C3*exp(_c[4]*t)*_C7*_C5*...
u(r,theta,z,t) = exp(_c[3]^(1/2)*z)*exp(_c[2]^(1/2)*theta)*_C3*exp(_c[4]*t)*_C7*_C5*_C1*BesselJ((-_c[2])^(1/2),(_c[3]-_c[4])^(1/2)*r)+exp(_c[3]^(1/2)*z)*exp(_c[2]^(1/2)*theta)*_C3*exp(_c[4]*t)*_C7*_C5*...
u(r,theta,z,t) = exp(_c[3]^(1/2)*z)*exp(_c[2]^(1/2)*theta)*_C3*exp(_c[4]*t)*_C7*_C5*_C1*BesselJ((-_c[2])^(1/2),(_c[3]-_c[4])^(1/2)*r)+exp(_c[3]^(1/2)*z)*exp(_c[2]^(1/2)*theta)*_C3*exp(_c[4]*t)*_C7*_C5*...
u(r,theta,z,t) = exp(_c[3]^(1/2)*z)*exp(_c[2]^(1/2)*theta)*_C3*exp(_c[4]*t)*_C7*_C5*_C1*BesselJ((-_c[2])^(1/2),(_c[3]-_c[4])^(1/2)*r)+exp(_c[3]^(1/2)*z)*exp(_c[2]^(1/2)*theta)*_C3*exp(_c[4]*t)*_C7*_C5*...

>   

Multivariate Integral Calculus

The VectorCalculus  package greatly simplifies integration of both real- and vector-valued functions over geometric regions in R^2  and R^3 , such as rectangles, parallelepipeds, circles, ellipses, spheres, triangles, tetrahedra, and regions bounded by curves and surfaces.  

Definite Integrals of Multivariate Functions over Regions and Solids

VectorCalculus  overloads Maple's int  routine to include options for specifying the geometry of the region of integration.

Compute Int(sin(x)*cos(y)*tan(z),V)  over the region x = 0 .. Pi, y = 0 .. Pi/2, z = 0 .. Pi/4

>    int( sin(x)*cos(y)*tan(z), [x, y, z] = Parallelepiped( 0..Pi, 0..Pi/2, 0..Pi/4 ) );

ln(2)

Compute Int(exp(-x^2-y^2),A)  over the circle with center `<,>`(0,0)  and radius r , and Int(exp(-x^2-y^2-z^2),V)  over all of R^3

>    int( exp(-x^2-y^2), [x, y] = Circle( <0, 0>, r ) ),
int( exp(-x^2-y^2-z^2), [x, y, z] = Sphere( <0, 0, 0>, infinity ) );

-Pi*exp(-r^2)+Pi, Pi^(3/2)

Find the volume of the solid bounded by the paraboloid z = x^2+y^2 , the plane x = y , and the surface y = x^2 .  We use the Region  option

>    int( 1, [x,y,z] = Region( 0..1, x^2..x, 0..x^2+y^2 ) );

3/35

Compute Int(x+y+z,S)  where S  is the tetrahedron with vertices `<,>`(0,0, 0), `<,>`(1,0, 0), `<,>`(0,1, 0), `<,>`(0,0, 1)

>    int( x+y+z, [x,y,z] = Tetrahedron( <0,0,0>, <1,0,0>, <0,1,0>, <0,0,1> ) );

1/8

Find the area of the section of the ellipse x^2/4+y^2/9 = 1  from the angle alpha  through the angle beta

>    int( x, [x,y] = Sector( Ellipse( x^2/4 + y^2/9 - 1 ), alpha, beta ) );

40*(1/(5*cos(beta)^2+4))^(3/2)*sin(beta)*cos(beta)^2+32*(1/(5*cos(beta)^2+4))^(3/2)*sin(beta)-40*(1/(5*cos(alpha)^2+4))^(3/2)*sin(alpha)*cos(alpha)^2-32*(1/(5*cos(alpha)^2+4))^(3/2)*sin(alpha)
40*(1/(5*cos(beta)^2+4))^(3/2)*sin(beta)*cos(beta)^2+32*(1/(5*cos(beta)^2+4))^(3/2)*sin(beta)-40*(1/(5*cos(alpha)^2+4))^(3/2)*sin(alpha)*cos(alpha)^2-32*(1/(5*cos(alpha)^2+4))^(3/2)*sin(alpha)

>   

Line and Surface Integrals of Vector-Valued Functions

VectorCalculus  has top-level functions for path, line and surface integrals.  You can specify the paths of integration with options such as LineSegments , Path , and Arc .

>    SetCoordinates( cartesian[x,y] );

cartesian[x,y]

Compute the line integral Int(F(x,y),S)  for F(x,y) = `<,>`(-x-y,x-y)  over the unit square, demonstrating that F  is a nonconservative vector field.

>    LineInt( VectorField( <-x-y, x-y> ), LineSegments( <0,0>, <1,0>, <1,1>, <0,1>, <0,0> ) );

2

Compute the line integral Int(F(x,y),S)  for F(x,y) = `<,>`(x^2,y^2)  over the path `<,>`(t,t^2)  as t  ranges from 0 to 2 .

>    LineInt( VectorField( <x^2, y^2> ), Path( <t, t^2>, t=0..2 ) );

24

Compute the line integral Int(F(x,y),S)  for F(x,y) = `<,>`(y,-x)  over the ellipse x^2/4+y^2/9 = 1   in the first quadrant .

>    LineInt( VectorField( <y,-x> ), Arc( Ellipse( x^2/4+y^2/9-1 ), 0, Pi/2 ) );

-3*Pi

>    restart: with(VectorCalculus):

Warning, the assigned names <,> and <|> now have a global binding

Warning, these protected names have been redefined and unprotected: *, +, ., Vector, diff, int, limit, series

Compute the surface area of a sphere

>    SurfaceInt( 1, [x,y,z] = Surface( <r,s,t>, s=0..Pi, t=0..2*Pi, coords=spherical ) )
assuming r>0;

4*Pi*r^2

Compute the surface area of the plane z = 4-2*x-y  over the triangle in the xy-plane that has vertices `<,>`(0,0), `<,>`(1,0), `<,>`(0,1)

>    SurfaceInt( 1, [x,y,z] = Surface( <s,t,4-2*s-t>, [s,t] = Triangle(<0,0>,<1,0>,<0,1>) ) );

1/2*6^(1/2)

>   

Operations on Parametric Curves

The VectorCalculus  package provides tools to analyse the differential-geometric properties of curves, such as torsion, normal vectors, radius of curvature, TNB frames, and others.

Example: Draw the evolute of an ellipse

The evolute  of a curve is the locus of its radius of curvature.  To compute the evolute of an ellipse, we first express the ellipse as a parametric curve.

>    c := <3*cos(t)/2, sin(t)>;

c := Vector(%id = 21756244)

>    assume( t::real ): interface( showassumed=0 ):

Compute the normal vector of the ellipse as a function of t .

>    n := simplify( PrincipalNormal( c, t ) );

n := Vector(%id = 2886996)

Scale the normal vector to a unit vector by dividing it by its length.

>    n_unit := simplify( n / LinearAlgebra :- Norm( n, 2) );

n_unit := Vector(%id = 17615492)

Compute the radius of curvature of the ellipse

>    r := RadiusOfCurvature( c, t );

r := -1/12*(5*cos(t)^2-9)*(-5*cos(t)^2+9)^(1/2)

The evolute is then the normal vector scaled by the radius of curvature at each value of t .

>    evolute := c + r * n_unit;

evolute := Vector(%id = 19488428)

>    plot([[c[1], c[2], t=0..2*Pi], [evolute[1], evolute[2], t=0..2*Pi]], scaling=constrained);

[Maple Plot]

Extending the VectorCalculus Package

You can add your own coordinate systems to the database of predefined systems in VectorCalculus using the AddCoordinates function .

>    restart: with(VectorCalculus): interface(showassumed=0):

Warning, the assigned names <,> and <|> now have a global binding

Warning, these protected names have been redefined and unprotected: *, +, ., Vector, diff, int, limit, series

Let's define the polarbear coordinate system   as
x(R,delta) = R^2*cos(delta), y(R,delta) = R^2*sin(delta)

>    assume( R >= 0, delta >=0, delta <=2*Pi ):

>    AddCoordinates( 'polarbear'[R,delta], [R^2*cos(delta),R^2*sin(delta)] );

polarbear

Now we can perform any computation we wish in polarbear coordinates.

>    Gradient( R^2*delta^2 - R/4, polarbear[R,delta] );

Vector(%id = 19584372)