Jacobian

Examples of Jacobian matrices

1. Multiple integrals
   In single integrations, integration by substitution (u-substitution) is stated as
   

   ⌠ ⌡ b a  f(x) d x = ⌠ ⌡ β α  f(x(u)) d x d u d u.

   (8)

   In multiple integrations, the above formula can be extended to
   
   

   ⌠ ⌡ ⌠ ⌡ D  f(x, y) d x d y = ⌠ ⌡ ⌠ ⌡ D′  f(x(u), y(v)) ⎢ ⎢ ∂(x, y) ∂(u, v) ⎢ ⎢ d u dv,

   (9)

   where the determinant of the Jacobian matrix,
   

   | J | ≡ ⎢ ⎢ ∂(x, y) ∂(u, v) ⎢ ⎢ ,

   is called the Jacobian determinant or simply the Jacobian. The Jacobian is hence interpreted as the scaling factor of an area element from one coordinate system to another.
   (Proof) Consider the area spanned by the two base vectors, →AB and →AC, below:
   When the coordinates of point A are denoted as (x, y), the coordinates of the points B and C are expressed as (x + ∂x/∂u du, y + ∂y/∂u du) and (x + ∂x/∂v dv, y + ∂y/∂v dv), respectively, so the components of the vectors, →AB and →AC, are
   
   

   → AB   = ⎛ ⎝ ∂x ∂u du, ∂y ∂u du ⎞ ⎠ , → AC   = ⎛ ⎝ ∂x ∂v dv, ∂y ∂v dv ⎞ ⎠ .

   Therefore, the area of the parallelogram spanned by the two vectors is ¹
   
   

   dS = ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ∂x ∂u , ∂x ∂v ∂y ∂u , ∂y ∂v ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ du dv = ⎢ ⎢ ∂(x, y) ∂(u, v) ⎢ ⎢ du d v = | J | d u dv.

   Exercise: Find |J| for the polar coordinate system:
   

   x = r cos θ,     y = r sin θ.

   (Answer)
   
   

   | J | = ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ∂x ∂r , ∂x ∂θ ∂y ∂r , ∂y ∂θ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ = ⎢ ⎢ ⎢ ⎢ ⎛ ⎜ ⎜ ⎝ cos θ, −r sin θ sin θ, r cos θ ⎞ ⎟ ⎟ ⎠ ⎢ ⎢ ⎢ ⎢ = r cos2 θ+ r sin2 θ = r.

   Exercise: Evaluate
   

   I ≡ ⌠ ⌡ ∞ −∞  e−x2 dx.

   (Answer)
   

   I2 = ⌠ ⌡ ∞ −∞  e−x2 dx ⌠ ⌡ ∞ −∞  e−y2 dy = ⌠ ⌡ ∞ −∞  ⌠ ⌡ ∞ −∞  e−(x2 + y2) dx dy = ⌠ ⌡ 2π 0  ⌠ ⌡ ∞ 0  e−r2 r dr dθ = 2π ⌠ ⌡ ∞ 0  e−r2 r dr = 2π× 1 2 = π,

   so
   

   I = √   π   ,

   where r² = t was used in evaluating ∫₀^r2 r e^−r2 dr .
2. Newton-Raphson method (single variable)
   The formula for the Newton-Raphson method can be derived by retaining the first two terms of the Taylor series of f(x) as
   

   f(x) = f(xo) + f′(xo) (x − xo) + (higher order terms).

   (10)

   With f(x) = 0, the above equation can be solved for x as
   

   x = xo − f(xo) f′(xo) ,

   (11)

   which can be used as the iteration scheme:
   For example, √2 can be computed by solving f(x) ≡ x² − 2 = 0.
   

   f(x) = x2 −2,     f′(x) = 2 x.

   Start with x₁ = 2.0 (initial guess), then,
   

   x1 = 2, x2 = 2− f(2) f ′(2) = 2 − 2 4 = 1.5, x3 = 1.5 − f(1.5) f ′(1.5) = 1.5 − 0.25 3 = 1.41667, x4 = 1.4167− f(1.4167) f ′(1.4167) = … = 1.4142.

   The convergence is attained after only 4 iterations.
   More examples:
   1. Square root of a ( = √a)
      Let f(x) = x² − a. It follows
      

      xn+1 = xn − xn2 − a 2 xn = 1 2 ⎛ ⎝ xn + a xn ⎞ ⎠ . (12)

      Example (√3 = 1.732 …)
      

      x1 = 1,     x2 = 1 2 (1 + 3 1 ) = 2,    x3 = 1 2 (2 + 3 2 ) = 1.75,    x4 = 1 2 (1.75 + 3 1.75 ) = 1.732.
   2. Inverse of a ( = 1/a)
      Let f(x) = 1/x − a. It follows
      

      xn+1 = xn − 1 xn −a −1 xn2 = xn (2 − a xn). (13)

      Example (1/6 = 0.16667…)
      

      x1 = 0.2,     x2 = 0.2 (2 − 6×0.2)=0.16,    x3 = 0.16 (2 − 6×0.16) = 0.1664, x4 = 0.1664 (2 − 6×0.1664) = 0.166666.
   3. 1/√a
      Let f(x) = 1/x² − a, it follows
      
      

      xn+1 = xn − 1/xn2 − a − 2/xn3 = xn (3 − a xn2)/2. (14)
3. Newton-Raphson method (multiple variables)
   We want to solve a set of (non-linear) simultaneous equations in the form of
   

   f1(x1, x2, x3, …xn) = 0, f2(x1, x2, x3, …xn) = 0, … … fn(x1, x2, x3, …xn) = 0.

   By expanding each of the above by the Taylor series, one obtains
   

   f1(x)  ∼ f1 (xo) + ∂f1 ∂x1 ⎢ ⎢ x0  (x1 − x10) + ∂f1 ∂x2 ⎢ ⎢ x0  (x2 − x20) + …+ ∂f1 ∂xn ⎢ ⎢ x0  (xn − xn0), f2(x)  ∼ f2 (xo) + ∂f2 ∂x1 ⎢ ⎢ x0  (x1 − x10) + ∂f2 ∂x2 ⎢ ⎢ x0  (x2 − x20) + …+ ∂f2 ∂xn ⎢ ⎢ x0  (xn − xn0), … fn(x)  ∼ fn (xo) + ∂fn ∂x1 ⎢ ⎢ x0  (x1 − x10) + ∂fn ∂x2 ⎢ ⎢ x0  (x2 − x20) + …+ ∂fn ∂xn ⎢ ⎢ x0  (xn − xn0).

   If x satisfies
   

   f1(x) = 0,     f2(x) = 0,     …fn(x) = 0,

   the above equations can be written as
   

   ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 : 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ f1 f2 : fn ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ + ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ∂f1 ∂x1 , ∂f1 ∂x2 , … ∂f1 ∂xn ∂f2 ∂x1 , ∂f2 ∂x2 , … ∂f2 ∂xn … … … … ∂fn ∂x1 , ∂fn ∂x2 , … ∂fn ∂xn ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ x1 − x10 x2 − x20 : xn − xn0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,

   (15)

   or
   
   

   0 = f(xo) + J (x − xo),

   (16)

   where J is the Jacobian matrix defined as
   
   

   J ≡ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ∂f1 ∂x1 ∂f1 ∂x2 … ∂f1 ∂xn ∂f2 ∂x1 ∂f2 ∂x2 … ∂f2 ∂xn … … … … ∂fn ∂x1 ∂fn ∂x2 … ∂fn ∂xn ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ∂(f1, f2, …, fn) ∂(x1, x2, …, xn) . (17)

   Equation (16) can be solved for x as
   
   

   x = xo − J−1 f(xo),

   (18)

   where J⁻¹ is the inverse matrix of J.
   Example:
   Solve numerically for
   

   x3 + y2 = 1,     x y = 1 2 .

   (19)

   (Solution) Let
   

   f1 ≡ x3 + y2 − 1,     f2 ≡ x y − 1 2 .

   It follows
   

   J = ∂(f1, f2) ∂(x, y) = ⎛ ⎜ ⎜ ⎝ 3 x2, 2 y y, x ⎞ ⎟ ⎟ ⎠ ,

   and
   

   J−1 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ x 3 x3 − 2 y2 , − 2 y 3 x3 − 2 y2 − y 3 x3 − 2 y2 , 3x2 3 x3 − 2 y2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,

   so
   

   J−1 f = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ x 3 x3 − 2 y2 , − 2 y 3 x3 − 2 y2 − y 3 x3 − 2 y2 , 3x2 3 x3 − 2 y2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎝ x3 + y2 − 1 x y − 1/2 ⎞ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ x4−x (y2+1)+y 3 x3−2 y2 4 x3 y−3 x2−2 y3+2 y 6 x3−4 y2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .

   (20)

   Hence the iteration scheme is
   

   ⎛ ⎜ ⎜ ⎝ xn+1 yn+1 ⎞ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎝ xn yn ⎞ ⎟ ⎟ ⎠ − ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ xn4−xn (yn2+1)+yn 3 xn3−2 yn2 4 xn3 yn−3 xn2−2 yn3+2 yn 6 xn3−4 yn2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .

   (21)

   Sample C code:

   #include <stdio.h>
   #include <math.h>
   int main()
   {
   double x=1.0, y=1.0;
   int i;
   
   for (i=0;i<=10;i++)
     { x = x -(pow(x,4) + y - x*(1 + y*y))/(3*pow(x,3) - 2*y*y);
       y = y -(-3*x*x + 2*y + 4*pow(x,3)*y - 2*pow(y,3))/(6*pow(x,3) - 4*y*y);
     }
   
   printf("%f %f\n", x, y);
   return 0;
   }
   
   

   Sample Matlab/Octave code:

   x=1; y=1;
   for i=0:10
    eqs=[ x*x*x+y*y-1; x*y - 1/2];
    jacob=[3*x*x 2*y; y x];
    right=inv(jacob)*eqs;
    x=x-right(1);
    y=y-right(2);
   end;
   
   fprintf('%f %f\n', x, y);
   
   

   Sample Mathematica code:

   FindRoot[{x^3 + y^2 == 1, x y == 1/2}, {{x, 1}, {y, 1}}]
   
   

   Online C compiler
   Online Matlab/Octave
   Wolfram Alpha
   Starting (x, y) = (1.0, 1.0), the convergence was reached at (x, y) = (0.877275, 0.569947) after only 4 iterations. Note that this is just one of multiple roots. It is necessary to try different initial guess to obtain other roots.

Footnotes: