Help for a "Surface and Flux integrals" program

[salvomic]

hi,
I need some help to implement for the Prime a simple CAS programs to calculate Surface Integrals and Flux Integrals (Surface integrals with vector fields)

Surface integral:
\[ \int _\sigma f dS = \iint_{A}f(\sigma _1(u,v),\sigma _2(u,v),\sigma _3(u,v))\sqrt{I_1^2+I_2^2+I_3^2}dudv \]

Flux integral:
\[
\int _\sigma F\cdot \mathbf{n} \, dS = \iint_{\sigma }(f_1(\sigma _1(u, v),\sigma _2(u, v),\sigma _3(u, v))I_1)+f_2(\sigma _1(u, v),\sigma _2(u, v),\sigma _3(u, v))I_2+f_3(\sigma _1(u, v),\sigma _2(u, v),\sigma _3(u, v))I_3) dudv
\]
(where \( F = (f_1, f_2, f_3) \))

Where also:
\[
I_1(u,v) = det\begin{pmatrix}
\frac{\partial \sigma _2}{\partial u}(u,v) & \frac{\partial \sigma _2}{\partial v}(u,v)\\
\frac{\partial \sigma _3}{\partial u}(u,v) & \frac{\partial \sigma _3}{\partial v}(u,v)\\
\end{pmatrix}
\, ;
I_2(u,v) = det\begin{pmatrix}
\frac{\partial \sigma _3}{\partial u}(u,v) & \frac{\partial \sigma _3}{\partial v}(u,v)\\
\frac{\partial \sigma _1}{\partial u}(u,v) & \frac{\partial \sigma _1}{\partial v}(u,v)\\
\end{pmatrix}
\, ;
I_3(u,v) = det\begin{pmatrix}
\frac{\partial \sigma _1}{\partial u}(u,v) & \frac{\partial \sigma _1}{\partial v}(u,v)\\
\frac{\partial \sigma _2}{\partial u}(u,v) & \frac{\partial \sigma _2}{\partial v}(u,v)\\
\end{pmatrix}

\]

Which are parts of: \( I = (I_1, I_2, I_3) \)

I would like to start (or follow) from this program in the Prime Software Library that has a simple syntax (for linear and curvilinear integrals):
INPUT 4 parameters: 1. a function (scalar / vectorial), 2. parametrisation of a curve, 3. lower bound, 4. upper bound;
and a control for arguments (2 or 3) and for the case there is no input (and then the program show a little help)...

I'd think to extend these concepts (from (curvi)linear integrals to surface and flux), using a parametrisation of the surface like "u, v, σ(u, v)"...
\[ \sigma (u, v) : \left\{\begin{matrix}
x_{1}=\sigma_{1}(u, v)) \\
x_{2}=\sigma_{2}(u, v)) \\
x_{3}=\sigma_{3}(u, v))
\end{matrix}\right.
\,\, , (u,v) \in A \subset \mathbb {R}^2
\]
But I've no clear idea at the moment, ehm :-)

Thank you in advance!
Salvo

EDIT: edited typo in laTeX formulas (thanks Alex!)...

[Arno K]

Perhaps you can use this:

Code:

#cas
sfint(w):=
BEGIN

local f:=w[1];
local g:=w[2];
local uvals:=w[3];
local vvals:=w[4];
local gg;
IF TYPE(w)≠4 THEN 
PRINT("Help");
PRINT("Usage:sfint('f(x,y,z)',[φ1(u,v),φ2(u,v),φ3(u,v)],[ulow,uhigh],[vlow,vhigh])")
PRINT("f can be a vector, too, i.e.:['x*y*z','x+y','x^2+z^3']");
return 0;
END;
f:=subst(f,[x,y,z]=g);
gg:=transpose(grad(g,[u,v]));
gg:=cross(col(gg,1),col(gg,2));
IF type(f)==DOM_LIST THEN
f:=DOT(f,gg);
ELSE
f:=f*ABS(gg);
END;
f:=int(f,u,uvals[1],uvals[2])
  return int(f,v,vvals[1],vvals[2]);
END;
#end

I tried it with : sfint({'x*y*z',[u+v,u-v,2*u-3*v],[1,2],[3,4]})
and: sfint({'x*y*z',[u+v,u-v,2*u-3*v],[1,2],[3,4]})
both seem to provide correct results, must do that by hand to check.
If you need help you may enter for the param w what you want instead of a list.
Arno

[salvomic]

thank you Arno!
I'm trying, but perhaps I put wrong syntax, as I get always 0, i.e. with z=x^2+y^2, x=u, y=v, z=u^2+y^2, [0,1], [0,2π]...
See image...
Maybe it would be better to use less syntax not using ' ' and the list. As in the linear integral where I used "argv" and "argc" to use something like:
intcur(3*x-4*y+z,[SIN(t), COS(t), t],0,1) or intlin([3*x, 4*y, z], [1, t, t^2], 0, 1)...

Here it could be simply: sfint(x^2+y^2,[u,v,u^2+v^2],[0,1],[0,2π]). What do you think about?

Salvo

[Arno K]

That is weird, must first check it, as there are always 4 arguments the usage of no list may be better.Without that list entry it seems difficult to force an error that then writes the help text, so we need an iferr around the whole computation.
Why now the error condition always is triggered, hm. Replace in IF: TYPE(w)≠4 by type(w)≠ DOM_LIST then this is gone and your example provides about approx 804.57.

[salvomic]

(11-03-2017 11:16 PM)Arno K Wrote:  That is weird, must first check it, as there are always 4 arguments the usage of no list may be better.Without that list entry it seems difficult to force an error that then writes the help text, so we need an iferr around the whole computation.
Why now the error condition always is triggered, hm. Replace in IF: TYPE(w)≠4 by type(w)≠ DOM_LIST then this is gone and your example provides about approx 804.57.

ok with DOM_LIST!
but my example is wrong, I must try another

It's correct this: sfint([x,-y,1],[u,v,u*v],[-1,1],[-1,1]) -> 4
Some examples I tried had integral bounds in polar form, so they wasn't right, sorry...

Ok, your program handles also Flux integral!

[Arno K]

Ok, I did it by hand, using your formula and entering the final integrals in the prime and got, guess, ..., 804.57

[Arno K]

IFERR seems not to work in cas-programs, commenting out lets the program run, with IFERR-block there is some syntax error somewhere.

[salvomic]

(11-03-2017 11:41 PM)Arno K Wrote:  Ok, I did it by hand, using your formula and entering the final integrals in the prime and got, guess, ..., 804.57

right!
please, read the post above, that I've edited now...
In my book some examples have polar form, so I was confused with the bounds...

That integral become from z=x^2+y^2, x=u, y=v, z=u^2+v^2 integrated in a circumference with center in 0,0 and radius 1...

Your code works!
Also is right:
sfint('z/(√(1+u^2+v^2))', [u,v,u*v],[0,1],[0,1]) -> 1/4

[Arno K]

(11-03-2017 11:23 PM)salvomic Wrote:  

(11-03-2017 11:16 PM)Arno K Wrote:  That is weird, must first check it, as there are always 4 arguments the usage of no list may be better.Without that list entry it seems difficult to force an error that then writes the help text, so we need an iferr around the whole computation.
Why now the error condition always is triggered, hm. Replace in IF: TYPE(w)≠4 by type(w)≠ DOM_LIST then this is gone and your example provides about approx 804.57.

ok with DOM_LIST!
but my example is wrong, I must try another

It's correct this: sfint([x,-y,1],[u,v,u*v],[-1,1],[-1,1]) -> 4
Some examples I tried had integral bounds in polar form, so they wasn't right, sorry...

Ok, your program handles also Flux integral!

Checked that by hand, 4 is ok, so my program works like it shall and the only remaining problem is catching an error to provide help, I will think about that, tomorrow as it is late now. But you can skip the '' around the entries, that leaves only the {} around that all and you can provide help...

[salvomic]

(11-04-2017 12:06 AM)Arno K Wrote:  Checked that by hand, 4 is ok, so my program works like it shall and the only remaining problem is catching an error to provide help, I will think about that, tomorrow as it is late now. But you can skip the '' around the entries, that leaves only the {} around that all and you can provide help...

I saw we can avoid also {}:
sfint(3*x+y,[u,v,u*v],[0,1],[0,2π]) works, like the other form (for flux) with a vector at first place. So the syntax is a bit easier.
For the control, you are right; however already in some cases the Prime complaint with a warning "no control is made for singular points..."

Salvo

EDIT:
also they are correct these examples (you can use also these ones for testing):
sfint(1,[r*COS(u), r*SIN(u), v],[0,2π],[0,L]) -> 2*L*π*|r|, that's the area of a cylinder with radius r and high L
and
sfint([y, -x, z],[r*SIN(u)*COS(v), r*SIN(u)*SIN(v), r*COS(u)],[0,π],[0,2π] ), the Flux of f(x)=yi-xj+zk exiting from a sphere radius 1 centered in the origin, with spheric coordinates in [0,π] x [0,2π] where σ(u,v) = r*sin(u)*cos(v) i + r*sin(u)*sin(v)j+r*cos(u)k
\( \int_{\Sigma }\mathbf{f}\cdot \mathbf{n}=\frac{4}{3}\pi r^3 \)

[salvomic]

I can't solve this problem; I'm following a book that perhaps report a wrong result: (π/6)*(5^(3/2)-1) -> 5.33041350027

Area of a paraboloid with cartesian equation z=x^2+y^2, integrated in a circumference with center in 0,0 and radius 1; we have x=u, y=v, z=u^2+v^2
We want calculate the area of the surface, so the function treated is constant 1.
The book suggest
\[ \int_{\sigma }1dS=\iint\sqrt{4u^2+4v^2+1}dudv=\int_{0}^{1}\int_{0}^{2\pi }\rho \sqrt{1+4\rho ^2}d \rho d\theta = \frac{\pi }{6} (5^\frac{3}{2}-1) \]
passing to the polar coordinates...
The Prime return another number for the integral (see attached image) and not \( \frac{\pi }{6} (5^\frac{3}{2}-1) \) that's about 5.33041350027.
Prime return 166.85...
So my previous example was so strange to test sfint()...
sfint(1,[u,v,u^2+v^2],[0,1],[0,1]) returns another number: 1.861564...

Salvo

[AlexFekken]

Hi Salvo (and Arno),

Nice to see some higher mathematics being implemented: I sincerely appreciate your work, so apologies if this seems like I am doing nothing more than nit-picking.

These are some typos that you might want to fix (and suggestions that you may or may not want to take up) if this is going to be the starting point for documentating your library:

1:

(11-03-2017 05:53 PM)salvomic Wrote:  \[ \int _\sigma dS = \iint_{A}f(\sigma _1(u,v),\sigma _2(u,v),\sigma _3(u,v))\sqrt{I_1^2+I_2^2+I_3^2}dudv \]

The f is missing on the left hand side:

\[ \int _\sigma f dS = \iint_{A}f(\sigma _1(u,v),\sigma _2(u,v),\sigma _3(u,v))\sqrt{I_1^2+I_2^2+I_3^2}dudv \]

And why not allow for vector functions f here as well? I realise that it is less common than a flux integral but the (component-wise) implementation would be trivial. Also, in my opinion, disitnguishing between 'regular' surface integrals and flux integrals more explicitly (e.g. by using a different function, like you did with your line and curvilinear integrals) might be a good thing.

2:

(11-03-2017 05:53 PM)salvomic Wrote:  \[
\int _\sigma F\cdot \mathbf{n} \, dS = \iint_{\sigma }(f_1(\sigma _1(u, v),\sigma _2(u, v),\sigma _3(u, v))I_1)+f_2(\sigma _1(u, v),\sigma _2(u, v),\sigma _3(u, v))I_2+f_3(\sigma _1(u, v),\sigma _2(u, v),\sigma _3(u, v))I_3) dudv
\]

Perhaps note that \[F = (f_1, f_2, f_3)\] for those less familiar with the maths.

3:

(11-03-2017 05:53 PM)salvomic Wrote:  \[ (u,v) \in A \supset \mathbb{R}^2 \]

The inclusion is the wrong way around:
\[ (u,v) \in A \subset \mathbb{R}^2 \]

You may also want to mention that A should be a rectangle, at least for the current implementation.

4:
Not really a typo, but for a better connect to the literature as well as the code it might be worth pointing out that \[I = (I_1, I_2, I_3)\] is the cross-product of the derivatives of \[ \sigma = (\sigma_1, \sigma_2, \sigma_3) \] w.r.t. u and v.

[salvomic]

thanks Alex,
I edited the first post with your corrections (they were good on the hook: I copied bad, sorry)...

I think we could separate the formulas for Flux and Surface or better specify them in the help.
Also I'm having trouble reproducing examples from books where, for facility purpose, they convert part of the double integrals into polar form: our program doesn't handle it well...

Also we could return a list with the I1,I2,I3 values (items of column vector of N(P), normal vector of the surface)...

However, Arno has already done 85% of works, already (and he, if wanted, could think to put a first version upon the Prime Software Library...) :-)

[AlexFekken]

(11-04-2017 09:14 AM)salvomic Wrote:  I can't solve this problem; I'm following a book that perhaps report a wrong result: (π/6)*(5^(3/2)-1) -> 5.33041350027

Hi Salvo,

You must paramtrise with polar coordinates from the start; the way you are doing it, your are not integrating over a circle in the (x, y)-space but a square (your image is missing but that is what I think you are doing)

E.g. parametrise the surface as

\[ (x=u \cos(v),\ y = u \sin(v),\ z= u^2)\]
where
\[(u,v) \in [0,1] \times [0,2 \pi] \]

and this will make

\[ I = (\cos(v), \sin(v), 2 u) \times (-u \sin(v), u \cos(v), 0) = (-2 u^2 \cos(v), - 2 u^2 \sin(v), u) \]

and

\[ |I| = \sqrt{4 u^4 + u^2} = u \sqrt{1 + 4 u^2} \]

etcetera...

[salvomic]

(11-04-2017 10:21 AM)AlexFekken Wrote:  Hi Salvo,

You must paramtrise with polar coordinates from the start; the way you are doing it, your are not integrating over a circle in the (x, y)-space but a square (your image is missing but that is what I think you are doing)

E.g. parametrise the surface as

\[ (x=u \cos(v),\ y = u \sin(v),\ z= u^2)\]
where
\[(u,v) \in [0,1] \times [0,2 \pi] \]

...

\[ |I| = \sqrt{4 u^4 + u^2} = u \sqrt{1 + 4 u^2} \]

etcetera...

It's true, Alex!
the correct input is, therefor:
sfint(1,[u*COS(v),u*SIN(v), u^2],[0,1],[0,2π]) that gives (1/6)*(√5*5*π-π) = 5.3304135, and simplifying we have the same result with (π/6)*(5^(3/2)-1) as in the book (that's correct!)

This fact induce me to think that we should find a way to help user with the changes of coordinates too...

EDIT: add:
But, if the result now correspond a that of the book, why in the Prime, doing this integral \( \int_{\sigma }1dS=\iint\sqrt{4u^2+4v^2+1}dudv=\int_{0}^{1}\int_{0}^{2\pi }\rho \sqrt{1+4\rho ^2}d \rho d\theta = \frac{\pi }{6} (5^\frac{3}{2}-1) \) we get a different result? Am I wrong with its input? or the book's example is not so clear...

Salvo

[AlexFekken]

(11-04-2017 10:08 AM)salvomic Wrote:  Also we could return a list with the I1,I2,I3 values (items of column vector of N(P), normal vector of the surface)...

Not sure what you expect from that but keep in mind that while I is a normal vector in the sense that it is orthogonal to the surface, it is not normalised to length 1 (unlike n). The relation between the differentials is

\[\mathbf{I} du dv = \mathbf{n} dS \]

and in this expression only n is normalised. This means that the length of I depends on the chosen parametrisation in terms of (u, v), e.g. you can scale (u,v) up or down and I will change accordingly.

In short: if you want to return a normal vector that is also normalised (length 1) then you must divide I by its l2 norm.

[salvomic]

(11-04-2017 10:40 AM)AlexFekken Wrote:  

(11-04-2017 10:08 AM)salvomic Wrote:  Also we could return a list with the I1,I2,I3 values (items of column vector of N(P), normal vector of the surface)...

Not sure what you expect from that but keep in mind that while I is a normal vector in the sense that it is orthogonal to the surface, it is not normalised to length 1 (unlike n). The relation between the differentials is

\[\mathbf{I} du dv = \mathbf{n} dS \]

and in this expression only n is normalised. This means that the length of I depends on the chosen parametrisation in terms of (u, v), e.g. you can scale (u,v) up or down and I will change accordingly.

In short: if you want to return a normal vector that is also normalised (length 1) then you must divide I by its l2 norm.

ok, I understand, it's so, indeed. Maybe return the I parts isn't so important if we don't give a way to handle the different parametrisation separately...

Please, Alex, read also my previous post, in the EDIT part, thank you.

[AlexFekken]

(11-04-2017 10:38 AM)salvomic Wrote:  This fact induce me to think that we should find a way to help user with the changes of coordinates too...

I disagree: the formulas take care of that automatically provided that you parametrise your surface (or curve, or manifold) correctly.

Hence my earlier remark that A needs to be a rectangle. So if the range of integration in your "natural" corrdinate space is not a rectangle then you have got your parametrisation wrong and need to fix that first before doing the integration.

You might want to introduce a separate function for transforming the differentials between one parametrisation and another, but I thnk the Prime already has an implementation of the required Jacobian that can be used...

[AlexFekken]

(11-04-2017 10:38 AM)salvomic Wrote:  EDIT: add:
But, if the result now correspond a that of the book, why in the Prime, doing this integral \( \int_{\sigma }1dS=\iint\sqrt{4u^2+4v^2+1}dudv=\int_{0}^{1}\int_{0}^{2\pi }\rho \sqrt{1+4\rho ^2}d \rho d\theta = \frac{\pi }{6} (5^\frac{3}{2}-1) \) we get a different result? Am I wrong with its input? or the book's example is not so clear...

Salvo

What exactly is the calculation that you carried out on the Prime that gave an incorrect result?

[Arno K]

If I use f(x,y)=x^2+y^2 and substite x=u*cos(v), y=u*sin(v),z=u^2 and integrate u from 0 to 1, v=0 to 2*PI I get (1/60)*(sqrt(5)*25*π+π), the same is provided by the prime, so some thing with that substitution must be wrong. Clearly my integral then has u^3*sqrt(...) which explains the difference, but where does it result from.