7.136 bug in Vector Valued Functions in Maple 7 and workaround (23.11.01)

7.136.1 PierLuigi Zezza
7.136.2 Robert Israel (26.11.01)
7.136.3 Wilhelm Werner (27.11.01)
7.136.4 Preben Alsholm (27.11.01)
7.136.5 Dr Francis J. Wright (28.11.01)
7.136.6 Koch-Beuttenmueller (29.11.01)
7.136.7 PierLuigi Zezza (29.11.01)
7.136.8 Robert Israel (3.12.01)

7.136.1 PierLuigi Zezza

In previous versions of Maple I was used to write vector valued functions in the following way

 > h:=vector(2): 
 > h[1]:=(x,y)->y-x^2: 
 > h[2]:=(x,y)->-(x)^2-(y+1)^2+1: 
 > h(x,y); 
 
                      [     2    2          2    ] 
                      [y - x , -x  - (y + 1)  + 1]
 

and it was very convenient because I could use single components of the function h or I could do things like

f:=(x,y,z)-> x^2+y^2-3*z^2; 
                                       2    2      2 
                    f := (x, y, z) -> x  + y  - 3 z 
 
 > Df:=vector(3,i->D[i](f)); 
 
     Df := [(x, y, z) -> 2 x, (x, y, z) -> 2 y, (x, y, z) -> -6 z]
 

With maple 7 it does not work any longer. If I try to use again the function h, I obtain

 > h(x,y); 
 
          [                 2          2                2    ] 
          [y(x, y) - x(x, y) , -x(x, y)  - (y(x, y) + 1)  + 1]
 

or in the second example

 > Df(x,y,z); 
 
                            [2 x, 2 y, -6 z] 
 
 > Df(1,2,3); 
 
              [2 x(1, 2, 3), 2 y(1, 2, 3), -6 z(1, 2, 3)]
 

I would like to understand what is changed and to find an alternative way of handling vector valued functions

It is corrected with Maple 8 (U. Klein)

7.136.2 Robert Israel (26.11.01)

It’s a very strange bug. Right after defining h I get:

> h(x,y); 
                               [     2     2          2    ] 
                               [x - y , - x  - (y + 1)  + 1]
 

So far it looks OK, but if you examine h:

> eval(h); 
                [               2              2          2    ] 
                [(x, y) -> y - x , (x, y) -> -x  - (y + 1)  + 1]
 

Somehow the definition of h has been changed so that it is the last value returned by h, not a vector of functions anymore. And then if you try h again, you get nonsense:

> h(s,t); 
                 [                 2           2                2    ] 
                 [x(s, t) - y(s, t) , - x(s, t)  - (y(s, t) + 1)  + 1]
 

It goes without saying that such a thing should never happen in Maple. I would suggest using lists, rather than vectors, of functions. Thus

> h:= [  (x,y) -> x - y^2, (x,y) -> -(x)^2-(y+1)^2+1 ]: 
 
> h(x,y); 
                               2    2          2 
                         [x - y , -x  - (y + 1)  + 1] 
 
> h(s,t); 
                               2    2          2 
                         [s - t , -s  - (t + 1)  + 1]
 

7.136.3 Wilhelm Werner (27.11.01)

Similar problem: I previously wrote e.g.

H:=map(unapply,grad(f(x,y),[x,y]),x,y);
 

(similar for jacobian, hessian etc.) in order to get the gradient as a vector valued function. In Maple 7 the situation is: it sometimes works, sometimes it does’nt, I have not been able to figure out what goes on.

7.136.4 Preben Alsholm (27.11.01)

I tried a similar example on a ’Vector’, a ’Matrix’, a ’vector’, and a ’matrix’ in R7. ’Vector’ returns an escaped local x at the definition. This must be a bug and not a feature.

’Matrix’, ’vector’, and ’matrix’ seem to work OK when mapping ’apply’ is done.

’vector’ and ’matrix’ seem to get redefined when direct application is attempted.

This is what I tried:

First a Vector:

> V:=Vector([s->s^2,s->cos(s+1)]); # A local x has escaped 
 
                                    [x] 
                               V := [ ] 
                                    [x] 
 
> S:=indets(V); 
 
                               S := {x} 
 
> op(S)-x; 
 
                                x - x 
 
> type(op(S),`local`); 
 
                                 true
 

Now a Matrix:

> M:=Matrix([[s->s^2,s->cos(s+1)],[sin,cos]]); 
 
                       [      2                   ] 
                  M := [s -> s     s -> cos(s + 1)] 
                       [                          ] 
                       [  sin            cos      ] 
 
> M(7); 
 
                   [      2                   ] 
                   [s -> s     s -> cos(s + 1)](7) 
                   [                          ] 
                   [  sin            cos      ]
 

The argument is not automatically mapped, but it can be done explicitly:

> map(apply,M,7); 
 
                          [  49      cos(8)] 
                          [                ] 
                          [sin(7)    cos(7)] 
 
> restart;
 

Now a vector:

> v:=vector([s->s^2,s->cos(s+1)]); 
 
                        [      2                 ] 
                   v := [s -> s , s -> cos(s + 1)]
 

Mapping seems to work just fine:

> map(apply,v,7); 
 
                             [49, cos(8)] 
 
> map(apply,v,t); 
 
                           [ 2            ] 
                           [t , cos(t + 1)] 
 
> map(apply,v,r); 
 
                           [ 2            ] 
                           [r , cos(r + 1)]
 

Direct application seems OK to start with, but is destructive:

> v(t); 
 
                           [ 2            ] 
                           [t , cos(t + 1)] 
 
> map(apply,v,7); 
 
                        [    2               ] 
                        [t(7) , cos(t + 1)(7)] 
 
> v(r); 
 
                        [    2               ] 
                        [t(r) , cos(t + 1)(r)] 
 
> restart;
 

Now a matrix:

> m:=matrix([[s->s^2,s->cos(s+1)],[sin,cos]]); 
 
                           [m[1, 1]    m[1, 2]] 
                      m := [                  ] 
                           [  sin        cos  ]
 

Again mapping is fine:

> map(apply,m,7); 
 
                          [  49      cos(8)] 
                          [                ] 
                          [sin(7)    cos(7)] 
 
> map(apply,m,t); 
 
                        [   2                ] 
                        [  t       cos(t + 1)] 
                        [                    ] 
                        [sin(t)      cos(t)  ] 
 
> map(apply,m,r); 
 
                        [   2                ] 
                        [  r       cos(r + 1)] 
                        [                    ] 
                        [sin(r)      cos(r)  ]
 

but direct application is destructive:

> m(t); 
 
                        [   2                ] 
                        [  t       cos(t + 1)] 
                        [                    ] 
                        [sin(t)      cos(t)  ] 
 
> m(r); 
 
                     [      2                   ] 
                     [  t(r)       cos(t + 1)(r)] 
                     [                          ] 
                     [sin(t)(r)      cos(t)(r)  ] 
 
> map(apply,m,9); 
 
                  [         2                      ] 
                  [  t(r)(9)       cos(t + 1)(r)(9)] 
                  [                                ] 
                  [sin(t)(r)(9)      cos(t)(r)(9)  ]
 

In line 17 of the code for Vector the local variable initializer is defined. It seems that if an ’eval’ is put around one of the x’es there, things work:

> restart; 
> arg:=[s->s^2,s->cos(s+1)]:
 

Here is a much shortened version of initializer:

> initializer := [seq(`if`(not type(x,{'list', 'Vector', 'Matrix', 'array', 
'Array'}),x,uha),x = arg)]; 
 
                        initializer := [x, x]
 

Notice the following:

> seq(x, x=arg); 
 
                                 x, x 
 
> map(eval,%); 
 
                                 x, x
 

It helps to replace the one x with eval(x):

> seq(eval(x), x=arg); 
 
                             2 
                       s -> s , s -> cos(s + 1) 
 
> initializer2 := [seq(`if`(not type(x,{'list', 'Vector', 'Matrix', 'array', 
'Array'}),eval(x),uha),x = arg)]; 
 
                                     2 
              initializer2 := [s -> s , s -> cos(s + 1)]
 

It don’t know if this change would have undesirable consequences.

7.136.5 Dr Francis J. Wright (28.11.01)

I suggest that you need to use literally "vector-valued functions" rather than "vectors of functions". The following seems to work in Maple 7:

> h := (x,y) -> vector([y-x^2, x^2-(y+1)^2+1]); 
 
                             [     2   2          2    ] 
              h := (x, y) -> [y - x , x  - (y + 1)  + 1] 
 
> h(x,y); 
 
                     [     2   2          2    ] 
                     [y - x , x  - (y + 1)  + 1] 
 
> h(1,2); 
 
                               [1, -7] 
 
> f := (x,y,z) -> x^2+y^2-3*z^2; 
 
                                      2    2      2 
                   f := (x, y, z) -> x  + y  - 3 z 
 
> Df := () -> vector([D[1](f)(args),D[2](f)(args),D[3](f)(args)]); 
 
      Df := () -> [D[1](f)(args), D[2](f)(args), D[3](f)(args)] 
 
> Df(x,y,z); 
 
                           [2 x, 2 y, -6 z] 
 
> Df(1,2,3); 
 
                             [2, 4, -18]
 

or ...

> Df := () -> vector([seq(D[i](f)(args),i=1..3)]); 
 
             Df := () -> [seq(D[i](f)(args), i = 1 .. 3)] 
 
> Df(x,y,z); 
 
                           [2 x, 2 y, -6 z] 
 
> Df(1,2,3); 
 
                             [2, 4, -18]
 

The second definition of Df is equivalent to the first but slightly more elegant. I couldn’t find a satisfactory way to use the vector syntax that uses a mapping to initialize the vector, but it would not be significantly more elegant than using seq.

7.136.6 Koch-Beuttenmueller (29.11.01)

Wilhelm Werner wrote: ...

Do you mean the following difference between Maple7 and Maple6?:

Maple7:

> restart: 
> with(linalg): 
Warning, the protected names norm and trace have been redefined and 
unprotected 
 
> H:=map(unapply,grad(f(x,y),[x,y]),x,y); 
 
    H := [(x, y) -> diff(f(x, y), x), (x, y) -> diff(f(x, y), y)] 
 
> H1:=map(unapply,grad(sin(x)*cos(y),[x,y]),x,y); 
 
      H1 := [(x, y) -> cos(x) cos(y), (x, y) -> -sin(x) sin(y)] 
 
> H1(2*a,b); 
 
                 [cos(2 a) cos(b), -sin(2 a) sin(b)]  # 1. call 
 
> H1(2*a,b); 
 
 [cos(2 a)(2 a, b) cos(b)(2 a, b), -sin(2 a)(2 a, b) sin(b)(2 a, b)] #recurrences ? 
 
> H1(2*a,b)[1]; 
 
           cos(2 a)(2 a, b)(2 a, b) cos(b)(2 a, b)(2 a, b) #recurrences ? 
 
> H(x,y); 
 
                       [d           d         ] 
                       [-- f(x, y), -- f(x, y)]   # 1. call 
                       [dx          dy        ] 
 
> H(x,y); 
 
               [/d         \        /d         \      ] 
               [|-- f(x, y)|(x, y), |-- f(x, y)|(x, y)]   #recurrences ? 
               [\dx        /        \dy        /      ] 
 
> H(x,y)[1]; 
 
                       /d         \ 
                       |-- f(x, y)|(x, y)(x, y) #recurrences ? 
                       \dx        /
 

Maple6:

> restart: 
> 
> with(linalg): 
Warning, the protected names norm and trace have been redefined and 
unprotected 
 
> H:=map(unapply,grad(f(x,y),[x,y]),x,y); 
 
    H := [(x, y) -> diff(f(x, y), x), (x, y) -> diff(f(x, y), y)] 
 
> H1:=map(unapply,grad(sin(x)*cos(y),[x,y]),x,y); 
 
      H1 := [(x, y) -> cos(x) cos(y), (x, y) -> -sin(x) sin(y)] 
 
> H1(2*a,b); 
 
                 [cos(2 a) cos(b), -sin(2 a) sin(b)] # 1. call 
 
> H1(2*a,b); 
 
                 [cos(2 a) cos(b), -sin(2 a) sin(b)] # you get always the same result 
 
> H1(2*a,b)[1]; 
 
                           cos(2 a) cos(b) # you get always the same result 
 
> H(x,y); 
 
                       [d           d         ] 
                       [-- f(x, y), -- f(x, y)] # 1. call 
                       [dx          dy        ] 
 
> H(x,y); 
 
                       [d           d         ] 
                       [-- f(x, y), -- f(x, y)] # you get always the same result 
                       [dx          dy        ] 
 
> H(x,y)[1]; 
                              d 
                              -- f(x, y) # you get always the same result 
                              dx
 

7.136.7 PierLuigi Zezza (29.11.01)

Thanks for all suggestions. My final goal was to define first and second derivative of multivariate functions as functions, what I got is the following (still missing something on the second derivative)

Using list seems to be the best choice

 > alpha:=[t-> sin(t^2),t -> t^2]; 
 
                                        2         2 
                    alpha := [t -> sin(t ), t -> t ]
 

so that

 > alpha(t); 
 
                                   2    2 
                             [sin(t ), t ] 
 
 > D(alpha); 
 
                                   2 
                      [t -> 2 cos(t ) t, t -> 2 t] 
 
 > D(alpha)(t); 
 
                                   2 
                           [2 cos(t ) t, 2 t] 
 
 > D(alpha)(3); 
 
                             [6 cos(9), 6]
 

Or, may be less elegant,

 > alpha:=map(unapply,vector([sin(t^2),t^2]),t); 
 
                             [          2         2] 
                    alpha := [t -> sin(t ), t -> t ] 
 
 > D(alpha); 
 
                      [            2             ] 
                      [t -> 2 cos(t ) t, t -> 2 t] 
 
 > D(alpha)(t); 
 
                           [       2        ] 
                           [2 cos(t ) t, 2 t] 
 
 > D(alpha)(3); 
 
                             [6 cos(9), 6]
 

moreover if

 > f:=(x,y,z)-> x^2*y-y^3*sin(z)+z^4; 
 
                                   2      3           4 
                f := (x, y, z) -> x  y - y  sin(z) + z
 

then you can define

 > Df:=[seq(D[i](f),i=1..3)]; 
 
                                            2      2 
   Df := [(x, y, z) -> 2 x y, (x, y, z) -> x  - 3 y  sin(z), 
 
                        3             3 
         (x, y, z) -> -y  cos(z) + 4 z ]
 

so that

 > Df(x,y,z); 
 
                       2      2           3             3 
              [2 x y, x  - 3 y  sin(z), -y  cos(z) + 4 z ] 
 
 > Df(1,2,3); 
 
                  [4, 1 - 12 sin(3), -8 cos(3) + 108]
 

or in general

 > D1:=proc(f) local x,n,k,var: n:=nops({op(1,eval(f))}); 
 > if n=1 then D(f) else 
 > [seq(D[h](f),h=1..n)] fi; end; 
 
   D1 := proc(f) 
local x, n, k, var; 
     n := nops({op(1, eval(f))}); 
     if n = 1 then D(f) else [seq(D[h](f), h = 1 .. n)] end if 
end proc
 

so that

 > D1(f); 
 
                                      2      2 
   [(x, y, z) -> 2 x y, (x, y, z) -> x  - 3 y  sin(z), 
 
                        3             3 
         (x, y, z) -> -y  cos(z) + 4 z ] 
 
 > D1(f)(2,3,a); 
 
                                                     3 
                 [12, 4 - 27 sin(a), -27 cos(a) + 4 a ]
 

For the second derivative the solution is not completely satisfactory. If you define

 > D2:=proc(f)local n: n:=nops({op(1,eval(f))}); if n=1 then (D@@2)(f) 
 > else 
 > map(unapply,convert([D[1](D1(f)),D[2](D1(f)),D[3](D1(f))](op(1,eval(f) 
 > )),matrix),op(1,eval(f))) end if end; 
 
   D2 := proc(f) 
local n; 
     n := nops({op(1, eval(f))}); 
     if n = 1 then (D@@2)(f) 
     else map(unapply, convert( 
         [D[1](D1(f)), D[2](D1(f)), D[3](D1(f))]( 
         op(1, eval(f))), matrix), op(1, eval(f))) 
     end if 
end proc 
then 
 
 > D2(f)(x,y,z); 
 
               [2 y        2 x                 0        ] 
               [                                        ] 
               [                             2          ] 
               [2 x    -6 y sin(z)       -3 y  cos(z)   ] 
               [                                        ] 
               [           2            3              2] 
               [ 0     -3 y  cos(z)    y  sin(z) + 12 z ]
 

In Maple 8 you get: Error, (in convert/matrix) expecting array, rtable or list (U. Klein)

 > D2(f)(1,a,3); 
 
                [2 a         2                 0       ] 
                [                                      ] 
                [                            2         ] 
                [ 2     -6 a sin(3)      -3 a  cos(3)  ] 
                [                                      ] 
                [           2            3             ] 
                [ 0     -3 a  cos(3)    a  sin(3) + 108] 
 
 > D2(sin)(x); 
 
                                -sin(x)
 

but

 > D2(f); 
 
                    [?[1, 1]    ?[1, 2]       0   ] 
                    [                             ] 
                    [?[2, 1]    ?[2, 2]    ?[2, 3]] 
                    [                             ] 
                    [   0       ?[3, 2]    ?[3, 3]]
 

I do not understand what’s happening so I cannot fix it. If you do not want a matrix as an output but only a list this works fine.

 > DD2:=proc(f)local n: n:=nops({op(1,eval(f))}); if n=1 then (D@@2)(f) 
 > else 
 > [D[1](D1(f)),D[2](D1(f)),D[3](D1(f))] end if end; 
 
   DD2 := proc(f) 
local n; 
     n := nops({op(1, eval(f))}); 
     if n = 1 then (D@@2)(f) 
     else [D[1](D1(f)), D[2](D1(f)), D[3](D1(f))] 
     end if 
end proc 
and 
 
 > DD2(f); 
 
   [[(x, y, z) -> 2 y, (x, y, z) -> 2 x, 0], [(x, y, z) -> 2 x, 
 
                                                    2 
         (x, y, z) -> -6 y sin(z), (x, y, z) -> -3 y  cos(z)], [ 
 
                             2                       3              2 
         0, (x, y, z) -> -3 y  cos(z), (x, y, z) -> y  sin(z) + 12 z ] 
 
         ] 
 
 > DD2(f)(x,y,z); 
 
                                          2 
   [[2 y, 2 x, 0], [2 x, -6 y sin(z), -3 y  cos(z)], 
 
                 2          3              2 
         [0, -3 y  cos(z), y  sin(z) + 12 z ]] 
 
 > DD2(f)(1,2,3); 
 
   [[4, 2, 0], [2, -12 sin(3), -12 cos(3)], 
 
         [0, -12 cos(3), 8 sin(3) + 108]]
 

7.136.8 Robert Israel (3.12.01)

| PierLuigi Zezza wrote: For the second derivative the solution ...

It looks like Maple refuses to display function definitions in a matrix. For example:

> matrix([[x -> x^2, 1],[2,3]]); 
 
                                 [?[1, 1]    1] 
                                 [            ] 
                                 [   2       3]
 

There’s nothing wrong with the matrix, it just won’t display properly. This has to do with last-name evaluation in the `print/array/matrix` function.

On the other hand, Maple will display function definitions in a Matrix. So you can say

> Matrix(%); 
                                 [      2     ] 
                                 [x -> x     1] 
                                 [            ] 
                                 [   2       3]