The Chain Rule

Here are two differentiable mappings 
  F : R^3 -> R^2
  G : R^2 -> R^2

F[{x_, y_, z_}] = {x^2 + y^2 + z^2, 1/(x - y)}
G[{s_, t_}] = {E^s Cos[t], E^s Sin[t]}

...and here is their composition.

G[F[{x, y, z}]]

I can make Mathematica compute the Jacobian of F for me.  Note that I need the ‘Transpose’ command, because unfortunately Mathematica looks at matrices in terms of rows instead of columns

JF[{x_, y_, z_}] = Transpose[{D[F[{x, y, z}], x], D[F[{x, y, z}], y], D[F[{x, y, z}], z]}]
MatrixForm[JF[{x, y, z}]]

I also make Mathematica compute the Jacobian of G for me.  Note that I’m taking advantage of a little Mathematica shorthand this time, and by enclosing the whole command in ‘MatrixForm’ I make the output appear in human readable form only.

MatrixForm[JG[{s_, t_}] = D[G[{s, t}], {{s, t}}]]

Now I make Mathematica compute the derivative of the composition G*F in two different ways--once directly and once using the Chain Rule

A = MatrixForm[D[G[F[{x, y, z}]], {{x, y, z}}]]

MatrixForm[JG[F[{x, y, z}]].JF[{x, y, z}]]

Are those the same?  Fortunately, Mathematica will check that for us, too.

Simplify[A - B]

0

One important use of the chain rule occurs when one “iterates” a function.  Consider the following instance of the so-called “Henon” map

H[{x_, y_}] = {y, x + y^2 - 2}

Since R^2 is both the domain and range of H, we can compose H with itself as often as we want:

H[H[{x, y}]]
H[H[H[{x, y}]]]
H[H[H[H[{x, y}]]]]
Nest[H, {x, y}, 10];

...but we really wouldn’t want to have to differentiate the results.  

Happily, the chain rule says we don’t have to.  
Suppose for instance that I want to know the derivative of iterates of H at {x,y}={1,1} in the direction {1,-1}.  
I work as follows.  First I compute the Jacobian matrix of H

MatrixForm[JH[{x_, y_}] = D[H[{x, y}], {{x, y}}]]

(({{0, 1}, {1, 2 y}}))

Next I initilize
	n keeps track of how many times I’ve composed H with itself
	p[n] keeps track of the image of {1,1} under H^n
	JH[n] keeps track of the Jacobian matrix of H^n
	h[n] keeps track of the directional derivative of H^n

n = 0
MatrixForm[p[n] = {1, 1}]
MatrixForm[JH[n] = {{1, 0}, {0, 1}}]
MatrixForm[h[n] = {1,-1}]

Now I begin to iterate.  The key point is that by the chain rule JH^{n+1}(x,y) = JH(H^n(x,y)) * JH^n(x,y).

I can then execute the following code block n times in order to compute the derivative of H^n at {1,1} in the direction {1,-1}

MatrixForm[  p[n+1]=H[p[n]]                    ]
MatrixForm[  JH[n+1]=JH[p[n]].JH[n]   ]
MatrixForm[  h[n+1]=JH[p[n]].h[n]       ]
n=n+1