first, second and third derivative of a vector function

Question

I can define the following vector function

$f(x,y,z) = [x^2 - 1, x^3 + y^2, z]$

in MatLab or in Maple.

I want to find (and evaluate) the first, second and third derivatives of it. How to find : $\nabla{f}$, $\nabla^2{f}$ and $\nabla^3{f}$? Here vector function $f$ is differentiated with respect to vector $(x,y,z)$. I can do the simple job of finding $\nabla{f}$ with Maple as follows :

$f(x,y,z) := [x^2 - 1, x^3 + y^2, z]$

$\nabla{f} := Jacobian (f,[x,y,z])$

But how to differentiate $\nabla{f}$?

Is there a function in MatLab (or in Maple) which may take the above function as an input and evaluate its derivatives for a given value of (x,y,z)?

Thank you.

Answer 1

This is a solution to your problem.

• syms x y z %you are declaring the variables which we'll be used as arguments
• F=[x.^2-1;x.^3+y.^2;z] %this is your vector function
• Fx=[diff(f(1),x);diff(f(2),x);diff(f(3),x)] %this is $\frac{dF}{dx}$
• Fy=[diff(f(1),y);diff(f(2),y);diff(f(3),y)] %this is $\frac{dF}{dy}$
• Fz=[diff(f(1),z);diff(f(2),z);diff(f(3),z)] %this is $\frac{dF}{dz}$

Answer 2

Not sure about Matlab, but in Maxima you can differentiate the whole vector as diff([x.^2-1;x.^3+y.^2;z], x) e.g.. If you build a function which from any $G$ yields [diff(G,x); diff(G,y); diff(G,z)], you should be able to get the $k+1$ tensors $\nabla^k \mathbf{F}$ directly.

Answer 3

Hopefully I have understood your correctly.

The procedure T below will produce an Array of expressions. And T can be called again (recursively) to get the higher derivatives.

The unapply command can create appliable procedures from the results, as one way to evaluate subsequently at given points in (x,y,z) space.

restart:

f:=[x^3-1,x^3*y^2-y^3,z^2*x*y^2]; # alternatively Array([x^2-1,x^3-y^3,z])

                      [ 3       3  2    3   2    2]
                 f := [x  - 1, x  y  - y , z  x y ]


T:=proc(A::{Array,list},L::list(name))
    Array(op(ArrayTools:-Dimensions(A)),1..nops(L),
          ()->diff(A[args[1..nargs-1]],L[args[-1]]));
end proc:


ans1 := T(f,[x,y,z]); # Jacobian

                     [    2                           ]
                     [ 3 x           0           0    ]
                     [                                ]
                     [   2  2     3        2          ]
             ans1 := [3 x  y   2 x  y - 3 y      0    ]
                     [                                ]
                     [  2  2        2                2]
                     [ z  y      2 z  x y     2 z x y ]


func1 := unapply(ans1,[x,y,z]):

func1(3,5,7);

                         [  27     0     0]
                         [                ]
                         [ 675   195     0]
                         [                ]
                         [1225  1470  1050]

ans2 := T(ans1,[x,y,z]):

ans2[..,..,1],ans2[..,..,2],ans2[..,..,3];

      [ 6 x      0       0   ]  [  0         0          0   ]  
      [                      ]  [                           ]  
      [     2     2          ]  [   2       3               ]  
      [6 x y   6 x  y    0   ], [6 x  y  2 x  - 6 y     0   ], 
      [                      ]  [                           ]  
      [           2         2]  [   2         2             ]  
      [  0     2 z  y  2 z y ]  [2 z  y    2 z  x    4 z x y]  

        [  0        0       0   ]
        [                       ]
        [  0        0       0   ]
        [                       ]
        [     2                2]
        [2 z y   4 z x y  2 x y ]


func2 := unapply(ans2,[x,y,z]):

val := func2(3,5,7):

val[..,..,1], val[..,..,2], val[..,..,3];

          [ 18    0    0]  [  0    0    0]  [  0    0    0]
          [             ]  [             ]  [             ]
          [450  270    0], [270   24    0], [  0    0    0]
          [             ]  [             ]  [             ]
          [  0  490  350]  [490  294  420]  [350  420  150]


ans3 := T(ans2,[x,y,z]):

ans3[..,..,3,2], map(diff,ans2[..,..,3],y);

            [  0      0      0  ]  [  0      0      0  ]
            [                   ]  [                   ]
            [  0      0      0  ], [  0      0      0  ]
            [                   ]  [                   ]
            [4 z y  4 z x  4 x y]  [4 z y  4 z x  4 x y]

As mentioned, T can be called again, on earlier results. And a recursive procedure could be written to fold that process up into a new mechanism to produce the higher derivative results directly. Or T could be modified as follows,

R:=proc(A::{Array,list},L::list(name),n::posint) local i;
      Array(op(ArrayTools:-Dimensions(A)),seq(1..nops(L),i=1..n),
            ()->diff(A[args[1..nargs-n]],L[[args[-n..-1]]][]));
   end proc:

R(f, [x,y,z], 1); # Jacobian

R(f, [x,y,z], 2):
%[..,..,1], %[..,..,2], %[..,..,3];
R(f, [x,y,z], 3):
%[..,..,3,2], map(diff,R(f,[x,y,z],2)[..,..,3],y);

And those results from R could also be passed to unapply.