Home › Blogs › scumath's blog

A crossprod problem

scumath's picture
Posted on 2008-02-20 06:21 By scumath (Maple Rating 2 88)
Categories:

I want to plot the normal vector of the Mobius strip represented by the parametric equations:
x=(2+u*cos(v/2))*cos(v):

y=(2+u*cos(v/2))*sin(v):

z=u*sin(v/2):

It is well-known that the normal vector is the crossproduct of [diff(x(u,v),u),diff(y(u,v),u),diff(z(u,v),u)] and [diff(x(u,v),v),diff(y(u,v),v),diff(z(u,v),v)].

 

So I want Maple to calculate the partial derivatives and the crossproduct for me.

with(plots):with(linalg):

a:=2:b:=0.5:

x:=(u,v)->(a+u*cos(b*v))*cos(v):

y:=(u,v)->(a+u*cos(b*v))*sin(v):

z:=(u,v)->u*sin(b*v):

N:=crossprod([diff(x(u,v),u),diff(y(u,v),u),diff(z(u,v),u)],[diff(x(u,v),v),diff(y(u,v),v),diff(z(u,v),v)]):

surface:=plot3d([x(u,v),y(u,v),z(u,v)],u=-1..1,v=0..2*Pi):

curve:=spacecurve([x(0,v),y(0,v),z(0,v)],v=0..2*Pi,

thickness=2,color=red):

M1:=(u,v)->N[1](u,v): M2:=(u,v)->N[2](u,v): M3:=(u,v)->N[3](u,v):
M1(u,v);M2(u,v):M3(u,v);

K:=40:

for i from 1 to K do ti:=i*4*Pi/K:

normal_vector[i]:=spacecurve([x(0,ti)+s*M1(0,ti),y(0,ti)+s*M2(0,ti),z(0,ti)+s*M3(0,ti)],s=0..1,color=blue,thickness=4):

od:

A:=display(seq(normal_vector[i],i=1..K),insequence=true):

display(A,surface,curve,scaling=constrained);

It can be seen that the crossproduct N is OK, but the normal vector does not appear.

 

Then I calculate the partial derivatives and the crossproduct by hand.  ( N=[M1,M2,M3] below is the normal vector )

with(plots):

a:=2:b:=0.5:

x:=(u,v)->(a+u*cos(b*v))*cos(v):

y:=(u,v)->(a+u*cos(b*v))*sin(v):

z:=(u,v)->u*sin(b*v):
surface:=plot3d([x(u,v),y(u,v),z(u,v)],u=-1..1,v=0..2*Pi):

curve:=spacecurve([x(0,v),y(0,v),z(0,v)],v=0..2*Pi,

thickness=2,color=red):

M1:=(u,v)->b*cos(b*v)^2*sin(v)*u-sin(b*v)*(-b*u*sin(b*v)*sin(v)+(a+u*cos(b*v))*cos(v)):

M2:=(u,v)->sin(b*v)*(-b*u*sin(b*v)*cos(v)-(a+u*cos(b*v))*sin(v)-b*cos(b*v)^2*u*cos(v)):

M3:=(u,v)->cos(b*v)*cos(v)*(-b*u*sin(b*v)*sin(v)+(a+u*cos(b*v))*cos(v))-cos(b*v)*sin(v)*(-b*u*sin(b*v)*cos(v)-(a+u*cos(b*v))*sin(v)):

K:=40:

for i from 1 to K do ti:=i*4*Pi/K:

normal_vector[i]:=spacecurve([x(0,ti)+s*M1(0,ti),y(0,ti)+s*M2(0,ti),z(0,ti)+s*M3(0,ti)],s=0..1,color=blue,thickness=4) od:

A:=display(seq(normal_vector[i],i=1..K),insequence=true):

display(A,surface,curve,scaling=constrained);

And the normal vector appears.

I don't know what is wrong with the first program and how to solve the problem.

 

Comments

John Fredsted's picture

Use unapply

Comment on 2008-02-20 10:12 from John Fredsted (Maple Rating 4 566)

It seems to work if you define M1, M2, and M3 as follows:

M1 := unapply(N[1],(u,v));
M2 := unapply(N[2],(u,v));
M3 := unapply(N[3],(u,v));
scumath's picture

Thanks

Comment on 2008-02-20 11:53 from scumath (Maple Rating 2 88)

Yes, the new definition of M1, M2 and M3 works perfectly well.
Thank you very much.
But I don't quite understand the meaning of unapply.

with(plots): with(linalg):
a:=2:b:=0.5:
x:=(u,v)->(a+u*cos(b*v))*cos(v):
y:=(u,v)->(a+u*cos(b*v))*sin(v):
z:=(u,v)->u*sin(b*v):
N:=crossprod([diff(x(u,v),u),diff(y(u,v),u),diff(z(u,v),u)],[diff(x(u,v),v),diff(y(u,v),v),diff(z(u,v),v)]):
surface:=plot3d([x(u,v),y(u,v),z(u,v)],u=-1..1,v=0..2*Pi):
curve:=spacecurve([x(0,v),y(0,v),z(0,v)],v=0..2*Pi,
thickness=2,color=red):
M1 := unapply(N[1],(u,v)):
M2 := unapply(N[2],(u,v)):
M3 := unapply(N[3],(u,v)):
K:=80:
for i from 1 to K do ti:=i*4*Pi/K:
normal_vector[i]:=spacecurve([x(0,ti)+s*M1(0,ti),y(0,ti)+s*M2(0,ti),z(0,ti)+s*M3(0,ti)],s=0..1,color=blue,thickness=4) od:
A:=display(seq(normal_vector[i],i=1..K),insequence=true):
display(A,surface,curve,scaling=constrained);
 

John Fredsted's picture

An explanation, hopefully

Comment on 2008-02-20 13:41 from John Fredsted (Maple Rating 4 566)

It is one of those functions that I am not yet completely sure of myself how or when to use. I am quite certain that others here at mapleprimes can give far better explanations.

But after having identified in your code what seemed to be the problem area, I knew from experience that unapply might solve the problem. According to its help page unapply turns an expression into a functional operator. An example, closely related to your code:

N := u + v:
M_old := (u,v) -> N;
M_new := unapply(N,u,v);
M_old(s,t);
M_new(s,t);
                      M_old := (u, v) -> N
                    M_new := (u, v) -> u + v
                             u + v
                             s + t

Note that the parantheses used in my previous post are superfluous.

acer's picture

if we are lucky

Comment on 2008-02-20 14:01 from acer (Maple Rating 4 1410)

If we're lucky, someone will explain its relationship to the lambda calculus in detail.

Jacques has written in the past that, "Maple's unapply is the same as Church's lambda abstraction operator."

The help-page ?unapply says,

- The unapply command implements the lambda-expressions of lambda calculus.

For reference see ``An Implementation of Operators for Symbolic Algebra
Systems'' by G.H. Gonnet, SYMSAC July 1986.

What I wonder is whether the help ever claimed that, "The scoping behaviour of unbound names is not the same in the lambda calculus," and if so how that might be true.

acer

JacquesC's picture

Lambda Calculus

Comment on 2008-02-20 15:37 from JacquesC (Maple Rating 4 1699)

Usually the lambda calculus does not treat open terms, ie terms with unbound terms in them.  Only closed terms are given a semantic, usually.

Unbound names occuring in functions are usually treated as global.  [The only exception being when one explicitly inserts an escaped local from another scope into a function].  Since the usual lambda calculus does not possess mutable variables, never mind global 'variables', there is no need for trying to figure out what they might mean.

scumath's picture

unapply: a hard term

Comment on 2008-02-20 22:19 from scumath (Maple Rating 2 88)

It seems that unapply is too hard for an ordinary Maple user like me. But I wiil learn to try it whenever necessary.

scumath's picture

Good example

Comment on 2008-02-20 22:23 from scumath (Maple Rating 2 88)

Your example lets me better understand the meaning of unapply. Thank you.

scumath's picture

A comparison

Comment on 2008-02-20 22:41 from scumath (Maple Rating 2 88)

I use the following simpler example to show the difference between the two definitions of M1, M2 and M3. It seems that my original definition produces expressions of Mi, while the definition using unapply produces operators of Mi.

with(linalg):
x:=(u,v)->cos(v)*u:
y:=(u,v)->sin(v)*u:
z:=(u,v)->u+v:
N:=crossprod([diff(x(u,v),u),diff(y(u,v),u),diff(z(u,v),u)],[diff(x(u,v),v),diff(y(u,v),v),diff(z(u,v),v)]):
M1:=(u,v)->N[1](u,v);
M2:=(u,v)->N[2](u,v);
M3:=(u,v)->N[3](u,v);
M11 := unapply(N[1],u,v);
M22:= unapply(N[2],u,v);
M33 := unapply(N[3],u,v);
-------------------------------------------------------------------

                       M1 := sin(v) - cos(v) u

                       M2 := -sin(v) u - cos(v)

                                 2           2
                     M3 := cos(v)  u + sin(v)  u

                  M11 := (u, v) -> sin(v) - cos(v) u

                 M22 := (u, v) -> -sin(v) u - cos(v)

                                       2           2
                M33 := (u, v) -> cos(v)  u + sin(v)  u

>
 

John Fredsted's picture

A little problem, though

Comment on 2008-02-21 05:01 from John Fredsted (Maple Rating 4 566)

As in your original code you define M1, M2, and M3 as

M1 := (u,v) -> N[1](u,v);
M2 := (u,v) -> N[2](u,v);
M3 := (u,v) -> N[3](u,v);

which is problematic. I did not comment on it then, though, because I thought it would only obscure things. But inserting values (u,v) = (s,t), say, brings out the problem:

M1(s,t);
M2(s,t);
M3(s,t);
              sin(v)(s, t) - cos(v)(s, t) u(s, t)
              -sin(v)(s, t) u(s, t) - cos(v)(s, t)
                     2                       2        
         cos(v)(s, t)  u(s, t) + sin(v)(s, t)  u(s, t)

One might think (at least, that was what I did yesterday) that the following might work:

M1 := (u,v) -> N[1];
M2 := (u,v) -> N[2];
M3 := (u,v) -> N[3];

But it does not as the following clearly shows (no proper dependence on s and t):

M1(s,t);
M2(s,t);
M3(s,t);
                       sin(v) - cos(v) u
                       -sin(v) u - cos(v)
                           2           2  
                     cos(v)  u + sin(v)  u

This is usually where the unapply business becomes murky to me: these latter definitions of M1, M2, and M3 seem like functional operators (which is what unapply produces), so why do they not work?

Robert Israel's picture

scope of variables

Comment on 2008-02-21 13:33 from Robert Israel (Maple Rating 4 1555)

The value of N[1] is an expression containing the global variables u and v.

In the function definition

M1:= (u,v) -> N[1]

the u and v are formal parameters.  These have nothing to do with the global variables of the same name. 

When you call M1(s,t), every occurrence of u and v in the body of the function definition is replaced by s and t respectively.  But the body of the function definition does not actually contain u or v, it is only N[1].  This will have u and v when evaluated, but that evaluation hasn't taken place yet.  So there is nothing to replace, and s and t don't appear in the final result.

This is actually one of the main reasons for using unapply: 

M1 := unapply(N[1], (u,v));

will produce a function definition in which the current value of N[1], not just the name N[1], appears, and the u and v there will be the formal parameters u and v. 

acer's picture

mileage

Comment on 2008-02-21 14:36 from acer (Maple Rating 4 1410)

I get a lot of mileage out of trying to adhere to a simple principal of using function parameter names which do not coincide with global names in my expressions. It helps me keep it straight.

For example, you probably wouldn't have been puzzled at all if you'd written it as,

N1 := u+v;
M1 := (a,b) -> N1;

The unapply function is nice. Without it, we might be writing things like,

N1 := u+v;
M1 := (a,b)->subs({u=a,v=b},N1);
M1(s,t);

..and then we could get confused by this,

N1 := u+v;
M1 := (u,v)->subs({u=u,v=v},N1);
M1(s,t);

and resort to this,

N1 := u+v;
M1 := (u,v)->subs({:-u=u,:-v=v},N1);
M1(s,t);

And then we might realize that the above wouldn't work all inside a procedure in which u and v were locals. In desperation,

p:=proc() local u,v,N1,M1;
  N1 := u+v;
  M1 := subs({U=u,V=v},(u,v)->subs({U=u,V=v},N1));
  M1(s,t);
end proc:
p();

And then we'd go crazy.

acer

John Fredsted's picture

Scoping

Comment on 2008-02-22 04:49 from John Fredsted (Maple Rating 4 566)

Thanks, Israel and Acer.

It is certainly not the first time that I have had difficulties with these scoping issues, and, I am afraid, it will probably not be the last time. However, to minimize the risk of the latter I will bookmark your posts.

Comment viewing options

Select your preferred way to display the comments and click "Save settings" to activate your changes.
}