#
# 10.1   Vectors
# 10.1.1   Vector operations
#
#
#
# MAPLE SESSION 1
#
>    u := vector([1,-4,5]);
>    v := vector([2,3,7]);
>    u + v;
>    evalm(u + v);
>    evalm(u - v);
>    evalm(5*u - 3*u);
#
# 10.1.2   Length, dot product and cross product
#
#
#
# MAPLE SESSION 2
#
>    with(linalg): 
>    u := vector([x,y,z]);
>    v := vector([2,-5,6]);
>    norm(u);
>    norm(u,2);
>    norm(v,2);
#
#
# MAPLE SESSION 3
#
>    with(linalg): 
>    U := vector([u[1],u[2],v[2]]);
>    V := vector([v[1],v[2],v[2]]);
>    dotprod(U,V);
>    a := vector([1,2,3]);
>    b := vector([-3,5,7]);
>    dotprod(a,b);
#
#
# MAPLE SESSION 4
#
>    with(linalg): 
>    u := vector([2,1,2]); 
>    v := vector([1,1,0]); 
>    angle(u,v); 
>    simplify(%);
#
#
# MAPLE SESSION 5
#
>    with(linalg): 
>    u := vector([1,2,3]);
>    v := vector([5,-2,1]);
>    w := crossprod(u,v);
#
#
# MAPLE SESSION 6
#
>    dotprod(u,w); 
>    dotprod(v,w);
#
# 10.1.3   Plotting vectors
#
#
#
# MAPLE SESSION 7
#
>    with(plottools): 
>    with(plots): 
>    v := vector([2,3]);
>    vec := arrow([0,0],v,.1,.2,.2,color=red): 
>    display(vec):
#
#
# MAPLE SESSION 8
#
>    with(plots): 
>    with(plottools): 
>    with(linalg): 
>    u := normalize(vector([1,2,3])); 
>    v := normalize(vector([5,-2,1])); 
>    w := crossprod(u,v); 
>    uvec:=arrow([0,0,0],u,w,.1,0.2,.1,color=red): 
>    vvec:=arrow([0,0,0],v,w,.1,.2,.1,color=blue): 
>    wvec:=arrow([0,0,0],w,u,.1,.2,.1,color=green): 
>    utext:=textplot3d([u[1],u[2],u[3]," u "],color=black): 
>    vtext:=textplot3d([v[1],v[2],v[3]," v "],align=LEFT,
>  color=black): 
>    wtext:=textplot3d([w[1],w[2],w[3]," u x v "],align=LEFT,
>  color=black): 
>    c := sphere([0,0,0], 0.1,color=black): 
>    display(uvec,utext,vvec,vtext,wvec,wtext,c,
>  scaling=constrained,axes=boxed,orientation=[25,60]);
#
# 10.2   Lines and planes} 
# 10.2.1   Lines} 
#
#
#
# MAPLE SESSION 9
#
>    with(geom3d): 
>    point(P,1,2,3);
>    point(Q,4,-7,2);
#
#
# MAPLE SESSION 10
#
>    line(l,[P,Q]);
#
#
# MAPLE SESSION 11
#
>    Equation(l,t);
#
#
# MAPLE SESSION 12
#
>    P := vector([1,2,3]); 
>    Q := vector([4,-7,2]); 
>    PQ := evalm(Q - P);
#
#
# MAPLE SESSION 13
#
>    with(linalg): 
>    Q := vector([1,2,3]); 
>    P := vector([0,1,2]); 
>    L := vector([2,-3,5]); 
>    PQ := Q - P; 
>    LPQ := crossprod(L,PQ); 
>    dLPQ := norm(LPQ,2); 
>    dL := norm(L,2); 
>    dist := dLPQ/dL;
#
# 10.2.2   Planes
#
#
#
# MAPLE SESSION 14
#
>    with(geom3d): 
>    point(P0,1,1,1): 
>    point(P1,2,1,3): 
>    point(P2,3,2,1):
#
#
# MAPLE SESSION 15
#
>    plane(P,[P0,P1,P2],[x,y,z]): 
>    Equation(P);
#
#
# MAPLE SESSION 16
#
>    with(linalg): 
>    P := vector([0,0,3]); 
>    Q := vector([1,2,3]); 
>    N := vector([-2,4,1]); 
>    PQ := Q - P; 
>    NPQ := dotprod(N,PQ); 
>    dist := abs(NPQ)/norm(N,2);
#
# 10.3   Vector-valued functions
#
#
#
# MAPLE SESSION 17
#
>    with(linalg): 
>    F := vector([1, t, sqrt(t)]);
>    G := vector([sin(t), cos(t), t]);
>    crossprod(F,G);
#
#
# MAPLE SESSION 18
#
>    F + G; 
>    dotprod(F,G);
#
# 10.3.1   Differentiation and integration of vector functions
#
#
#
# MAPLE SESSION 19
#
>    F := vector([f[1](t),f[2](t),f[3](t)]);
>    map(diff,F,t);
#
#
# MAPLE SESSION 20
#
>    r := vector([2*cos(t), 3*sin(t), t]);
>    v := map(diff,r,t);
>    a := map(diff,v,t);
#
#
# MAPLE SESSION 21
#
>    r := vector([t^2, ln(1+t), sqrt(1-t)]); 
>    map(int,r,t); 
>    map(int,r,t=0..1);
#
# 10.3.2   Space curves
#
#
#
# MAPLE SESSION 22
#
>    with(plots): 
>    spacecurve([cos(t),3*sin(t),t],t=0..4*Pi,color=black,
>   thickness=3, numpoints=200, axes=boxed, 
>   orientation=[30,65]);
#
#
# MAPLE SESSION 23
#
>    fx:=(rho,phi,theta)->rho*sin(phi)*cos(theta);
>    fy:=(rho,phi,theta)->rho*sin(phi)*sin(theta);
>    fz:=(rho,phi,theta)->rho*cos(phi);
#
#
# MAPLE SESSION 24
#
>  plot3d([fx(1,phi,theta),fy(1,phi,theta),fz(1,phi,theta)],
>  phi=-Pi..Pi,theta=0..2*Pi);
#
#
# MAPLE SESSION 25
#
>    with(plots): 
>    S := spacecurve([cos(t),sin(t),t],t=0..4*Pi,color=black,
>  numpoints=200,axes=boxed): 
>    A:= animate3d([cos(t)+fx(1/10,phi,theta),sin(t)+
>  fy(1/10,phi,theta), t+ fz(1/10,phi,theta)],phi=-Pi..Pi,
>  theta=0..2*Pi,t=0..4*Pi, frames=32, 
>  scaling=constrained):
>    display(S,A,scaling=constrained,orientation=[60,30]);
#
#
# MAPLE SESSION 26
#
>    with(linalg): 
>    r := vector([cos(t),sin(t),t]);
>    dr := map(diff,r,t);
>    nr := simplify(norm(dr,2));
>    L := int(nr,t=0..4*Pi);
#
# 10.3.2   Tangents and normals to curves
#
#
#
# MAPLE SESSION 27
#
>    r := t -> 6*t*i + 3* sqrt(2)*t^2*j + 2*t^3*k;
#
#
# MAPLE SESSION 28
#
>    read vecfuncs: 
>    r := t -> 6*t*i + 3* sqrt(2)*t^2*j + 2*t^3*k: 
>    vec := r2vec(r(t));
#
#
# MAPLE SESSION 29
#
>    vec2r(vec);
#
#
# MAPLE SESSION 30
#
>    normv(vec);
#
#
# MAPLE SESSION 31
#
>    with(linalg): 
>    norm(2,vec);
#
#
# MAPLE SESSION 32
#
>    rp := diff(r(t),t);
#
#
# MAPLE SESSION 33
#
>    rpv := r2vec(rp);
>    n := normv(rpv);
#
#
# MAPLE SESSION 34
#
>    read vecfuncs: 
>    T := rp/n;
>    T := rsimp(T);
#
#
# MAPLE SESSION 35
#
>    normv(r2vec(T));
#
#
# MAPLE SESSION 36
#
>    nt := diff(T,t);
>    nt := rsimp(nt);
>    nn := normv(r2vec(nt));
>    N := nt/nn;
#
#
# MAPLE SESSION 37
#
>    normv(r2vec(N));
>    with(linalg): 
>    dotprod(r2vec(T),r2vec(N),orthogonal);
>    normal(%);
#
# 10.3.3   Curvature
#
#
#
# MAPLE SESSION 38
#
>    kappa := normal(nn/n);
#
#
# MAPLE SESSION 39
#
>    read vecfuncs: 
>    r := t -> 6*t*i + 3* sqrt(2)*t^2*j + 2*t^3*k: 
>    v := r2vec(diff(r(t),t)); 
>    a := r2vec(diff(r(t),t,t)); 
>    with(linalg): 
>    crossprod(v,a); 
>    normv(%)/normv(v)^3; 
>    normal(%);
#
# 10.4   The gradient and directional derivatives
#
#
#
# MAPLE SESSION 40
#
>    with(linalg): 
>    f := x^3 + sin(x +y*z^2);
>    grad(f, [x,y,z]);
#
#
# MAPLE SESSION 41
#
>    with(linalg): 
>    f := x*y+ 3*y*z^3;
>    v := vector([1,2,3]);
>    u := normalize(v);
>    grad(f,[x,y,z]);
>    g := subs(x=1,y=-1,z=2,%);
>    dotprod(g,u);
#
# 10.5   Extrema
# 10.5.1   Local extrema and saddle points
#
#
#
# MAPLE SESSION 42
#
>    f := (x,y) -> x^3 - 3*y*x + y^3;
>    criteqs := {diff(f(x,y),x)=0, diff(f(x,y),y)=0};
>    solve(criteqs, {x,y});
#
#
# MAPLE SESSION 43
#
>    f := (x,y) -> x^3 - 3*y*x + y^3;
>    with(linalg): 
>    h := hessian(f(x,y), [x,y]);
>    det(h);
>    des := unapply(%,x,y);
>    des(0,0);
>    des(1,1);
#
#
# MAPLE SESSION 44
#
>    f := (x,y) -> x^3 - 3*y*x + y^3; 
>    plot3d(f(x,y),x=-0.1..0.1,y=-0.1..0.1,axes=boxed,
>  style=patch,orientation=[20,70]);
#
#
# MAPLE SESSION 45
#
>    plot3d(f(x,y),x=0.9..1.1,y=0.9..1.1,axes=boxed,
>  style=patch,orientation=[20,70]);
#
# 10.5.2   Lagrange multipliers
#
#
#
# MAPLE SESSION 46
#
>    with(linalg): 
>    f := x^2 +y^2 + z^2;
>    g := x*y + x*z;
>    grad(f,[x,y,z])-lambda*grad(g,[x,y,z]);
>    evalm(%);
>    map(V->V=0,%): 
>    convert(%,set): 
>    EQNS := % union {g=4};
#
#
# MAPLE SESSION 47
#
>    \_EnvExplicit := true;
#
#
# MAPLE SESSION 48
#
>    SOL := solve(EQNS,x,y,z,lambda);
#
#
# MAPLE SESSION 49
#
>    subs(SOL[1],sqrt(f));
>    radsimp(%);
#
# 10.6   Multiple integrals
# 10.6.1   Double integrals
#
#
#
# MAPLE SESSION 50
#
>    I1 := Int(Int(sin(x^3),x=sqrt(y)..1),y=0..1);
>   value(%);
#
#
# MAPLE SESSION 51
#
>    I2:=Int(Int(sin(x^3),y=0..x^2),x=0..1);
>    value(%);
#
#
# MAPLE SESSION 52
#
>    evalf(I1); 
>    evalf(I2);
#
# 10.6.2   Triple integrals
#
#
#
# MAPLE SESSION 53
#
>    I1:=Int(z^2*sqrt(x^2+y^2),z=-sqrt(1-x^2-y^2)..
>  sqrt(1-x^2-y^2)): 
>    I2:=Int(I1,y=-sqrt(1-x^2)..sqrt(1-x^2)): 
>    I3:=Int(I2,x=-1..1);
>   value(%);
#
#
# MAPLE SESSION 54
#
>    f:=z^2*r;
>    I1:=Int(f*r,z=-sqrt(1-r^2)..sqrt(1-r^2)): 
>    I2:=Int(I1,r=0..1): 
>    I3:=Int(I2,theta=0..2*Pi);
>    value(I3);
#
#
# MAPLE SESSION 55
#
>    f :=z^2*r;
>    f := subs(r=rho*sin(phi),z=rho*cos(phi),f);
>    d := rho^2*sin(phi): 
>    I1:=Int(f*d,phi=0..Pi): 
>    I2:=Int(I1,rho=0..1): 
>    I3:=Int(I2,theta=0..2*Pi);
>    value(%);
#
# 10.6.3   The Jacobian
#
#
#
# MAPLE SESSION 56
#
>    r := rho*sin(phi): 
>    x := r*cos(theta): 
>    y := r*sin(theta): 
>    z := rho*cos(phi): 
>    f := vector([x,y,z]);
>    with(linalg): 
>    jacobian(f,[rho,phi,theta]);
>    det(%);
>    simplify(%);
#
# 10.7   Vector field
# 10.7.1   Plotting a vector field
# 10.7.2   Divergence and curl
#
#
#
# MAPLE SESSION 57
#
>    F := vector([x*y,y^2*z^2,x*y*z]);
>    with(linalg): 
>    diverge(F,[x,y,z]);
>    curl(F, [x,y,z]);
#
# 10.7.3   Potential functions
#
#
#
# MAPLE SESSION 58
#
>    with(linalg): 
>    F := vector([y*z+z^3-4*y^2, x*z-8*y*x, y*x+3*x*z^2]);
>    potential(F,[x,y,z],'Phi');
>    Phi;
#
#
# MAPLE SESSION 59
#
>    grad(Phi,[x,y,z])-F; 
>    evalm(%); 
>    simplify(%);
#
#
# MAPLE SESSION 60
#
>    with(linalg): 
>    B := vector([x-y,-y,-x]);
>    vecpotent(B,[x,y,z],'G');
>    evalm(G);
#
#
# MAPLE SESSION 61
#
>    curl(G,[x,y,z]);
#
# 10.8   Line integrals
#
#
#
# MAPLE SESSION 62
#
>    with(linalg): 
>    F := (x,y) -> vector([2*x*y,x-y]);
>    r := t -> vector([cos(t),sin(t)]);
>    dr := map(diff,r(t),t);
>    Int(dotprod(F(r(t)[1],r(t)[2]),dr,orthogonal),t=0..2*Pi);
>    value(%);
#
# 10.9   Green's Theorem
#
#
#
# MAPLE SESSION 63
#
>    M := 2*x*y;
>    N := x - y;
>    diff(N,x)-diff(M,y);
>    subs(x=cos(theta),y=sin(theta),%);
>    Int(Int(%*r,r=0..1),theta=0..2*Pi);
>    value(%);
#
# 10.10   Surface integrals
#
#
#
# MAPLE SESSION 64
#
>    x:=(3+cos(u))*cos(v); 
>    y:=(3+cos(u))*sin(v); 
>    z:=sin(u); 
>    plot3d([x,y,z],u=0..2*Pi,v=0..2*Pi,scaling=constrained,
>  orientation=[45,65]);
#
#
# MAPLE SESSION 65
#
>    with(linalg): 
>    x:=(R+r*cos(u))*cos(v);
>    y:=(R+r*cos(u))*sin(v);
>    z:=r*sin(u);
>    rv:=vector([x,y,z]);
>    Tu := map(diff,rv,u);
>    Tv := map(diff,rv,v);
>    cp := crossprod(Tu,Tv);
>    simplify(dotprod(cp,cp,orthogonal)): 
>    n := radsimp(sqrt(factor(%)));
>    int(int(n,u=0..2*Pi),v=0..2*Pi);
#
# 10.10.1   Flux of a vector field
#
#
#
# MAPLE SESSION 66
#
>    with(linalg): 
>    x:=sin(phi)*cos(theta);
>    y:=sin(phi)*sin(theta);
>    z:=cos(phi);
>    rv:=vector([x,y,z]);
>    Tphi := map(diff,rv,phi);
>    Ttheta := map(diff,rv,theta);
>    cp := simplify(crossprod(Tphi,Ttheta));
#
#
# MAPLE SESSION 67
#
>    F := (x,y,z) -> vector([y^4,z^4,x^4]);
>    dp := dotprod(F(x,y,z),cp,orthogonal);
>    int(int(dp,phi=0..Pi/2),theta=0..2*Pi);
#
# 10.10.2   Stoke's Theorem} 
#
#
#
# MAPLE SESSION 68
#
>    with(linalg): 
>    F := (x,y,z) -> vector([y^4,z^4,x^4]);
>    vecpotent(F(x,y,z),[x,y,z],G);
>    evalm(G);
#
#
# MAPLE SESSION 69
#
>    curl(G,[x,y,z]);
#
#
# MAPLE SESSION 70
#
>    r:=t->vector([cos(t),sin(t),0]); 
>    dr := map(diff,r(t),t);
>    rsubs := J -> subs({x=cos(t),y=sin(t),z=0},J);
>    GG := map(rsubs,G);
>    Gdr := dotprod(GG,dr,orthogonal);
>    int(Gdr,t=0..2*Pi);
#
# 10.10.3   The Divergence Theorem} 
#
#
#
# MAPLE SESSION 71
#
>    divF := 3*rho^2;
>    dV := rho^2*sin(phi);
>    int(int(int(divF*dV,theta=0..2*Pi),phi=0..Pi),rho=0..1);
#
#
# MAPLE SESSION 72
#
>    with(linalg): 
>    x:=sin(phi)*cos(theta): 
>    y:=sin(phi)*sin(theta): 
>    z:=cos(phi): 
>    rv:=vector([x,y,z]): 
>    Tphi := map(diff,rv,phi): 
>    Ttheta := map(diff,rv,theta): 
>    cp := simplify(crossprod(Tphi,Ttheta)): 
>    F := (x,y,z) -> vector([x^3,y^3,z^3]): 
>    dp := dotprod(F(x,y,z),cp,orthogonal): 
>    int(int(dp,phi=0..Pi),theta=0..2*Pi);