Question: How to represent functions when working with flux

Hi,

I'm trying to help the daughter of one of my friends with some flux related exercises.

We have calculated this one by hand and then I wanted to repeat it in Maple to help her see how to use Maple.

My problem is that I cannot figure out how to represent the functions in Maple, so that the calculations fall out similar to what we have been doing by hand. Yes, we could just use the built in Flux operation, but that defeats the purpose here.

Please bear with me - my calculus is rusty and I only have access to Maple on my son's machine. Plus I have never used Maple before... yes, I'm that old ;-)

We want to calculate the flux of the vector field V through the surface given by r:

with(Student[VectorCalculus]);

r:= (u,v) -> <u, v, 1-u^3>

V := VectorField(<x*z, x*cos(y), 3*x^2>

Then we can differentiate r wrt both u and v and calculate the normal vector for the surface:

Nf := ru &x rv

The next step is then to do the double integral of:

evalVF(V, r(u, v)) . Nf

This is - bar the evalVF call - almost identical to what you'd be doing by hand, so this is good.

Next is the flux out of the 3d object constrained by the square [0,1] and [-Pi/2, Pi/2] and the surface r.

The paramterisation of that is:

Omega := (u,v,w) -> <u, v, w*(1-u^3)>

with w in [0,1].

At this point one applies Gauss' Divergence Theorem which means that in order to find the flux of V out of Omega we "just" have to calculate the spatial integral of the divergence of V over the body described by Omega.

But that can be calculated by integrating V(Omega(u,v,w)) times the Jacobi function.

The Jacobi function is the determinant of the output of Maple's Jacobian procedure:

jacobi := Jacobian(Omega(u, v, w), [u, v, w])

In this case it is 1 - u^3.

Now comes the big problem: how do I get Div(V)(Omega(u,v,w)), ie, the divergence of V taken in the point Omega(u,v,w)?

If I do evalVF(V, Omega(u, v, w)) I get a column vector out that I cannot compute the divergence of using Divergence(%) since that only works on VectorFields.

And Divergence(V) is not a function so I cannot do

DivV := Divergence(V)

DivV(Omega(u,v,w))

And if I change Omega to return a value in R^3 instead of a vector I cannot use Jacobian, so that stop even faster.

Defining two versions of Omega just to make Maple happy seems wrong.

I just want to write the tripple integral of Div(V)(Omega(u,v,w))*jacobi or something close to that.

Please let me know if I should provide the workbook or in some other way make all my calculations available.

Thanks in advance and sorry if my description is lacking,

Torben