Another example of error analysis in numerical integration

Example

1. Left rectangular rule (5 points)
   
   

   I  ∼ 0.2 ( f(0) + f(0.2) + f(0.4) +f(0.6) + f(0.8)) = 2.36867.

   Way too underestimated (see the figure).
   
   
2. Mid point rule (5 points)
   
   

   I  ∼ 0.2 ( f(0.1) + f(0.3) + f(0.5) + f(0.7)+f(0.9)) = 3.1279.

   Much improved over the left rectangular rule.
   
   
3. Trapezoidal rule
   
   

   I  ∼  0.2 2 ( f(0) + 2 ×(f(0.2) + f(0.4) + f(0.6) + f(0.8)) + f(1.0) ) = 3.16867.

   The accuracy is the same as the midpoint rectangular rule.
   
   

Numerical differentiation

Some basics on infinite series

s=0;
for i=1:100 
 s=s+1.0/i;
end;

fprintf('sum=%f\n',s);

#include <stdio.h>
int main()
{
double s=0; int i;
for (i=1;i<100;i++) s=s + 1.0/i;
printf("sum=%f\n",s);
return 0;
}

The Riemann zeta (ζ) function (p-series)

Asymptotic notation and Landau-Bachmann magnitude (Big Oh)

Definition (Asymptotic notation)

Definition (Bachmann-Landau order of magnitude, aka Big Oh)

Examples

1. 
   

   3 x4 + x + 8 x2 + 2  ∼ 3 x2     as     x → ∞,
2. 
   

   3 x4 + x + 8 x2 + 2  ∼ 4     as     x → 0,
3. 
   

   3 x4 + x + 8 x2 + 2 = O(x2)     as     x → ∞,
4. 
   

   3 x4 + x + 8 x2 + 2 = O(x3)     as     x → ∞.

Convergence of infinite series

THEOREM (Comparison Test)

THEOREM (Comparison Test 2)

Counterexamples

THEOREM (Ratio Test)

Example

THEOREM (Alternating series) Leibniz's test (aka Dirichlet's test)

Example

Exercise of infinite series

1. 
   

   ∑ ⎛ ⎝ n2 + 2 n −1 n4 + 3 ⎞ ⎠ 3/2

   This series is convergent as a_n  ∼ 1/n³.
2. 
   

   ∑ 1 (n + 3)3/2

   This series converges as a_n  ∼ 1/n^1.5 (p > 1).
3. 
   

   ∑ n n2 + 3 ln n

   This series diverges as a_n  ∼ 1/n.
4. 
   

   ∑ ⎛ ⎝ 1 + 1 n2 ⎞ ⎠

   This series diverges as a_n ≠ 0 as n → ∞ ².
5. 
   

   ∑ ⎛ ⎝ cos n 2 n −1 ⎞ ⎠ 2

   This series converges as
   

   an ≤ 1 (2n − 1)2  ∼  1 4 n2 .
6. 
   

   ∑ (−1)n ln  ⎛ ⎝ 1 + 1 √n ⎞ ⎠

   This series converges per Leibniz's test, i.e ∑(−1)ⁿ is bounded and
   

   ln  ⎛ ⎝ 1 + 1 √n ⎞ ⎠ ↓ 0

   as n → ∞.
7. 
   

   ∑ 1 x2 + n2

   This series converges as a_n = O(1/n²).
8. 
   

   ∑ (−1)n 2n

   This series is convergent as |a_n+1/a_n| → 1/2 (ratio test).
9. 
   

   ∑ ln  ⎛ ⎝ n2 − 1 n2 + 1 ⎞ ⎠

   This series converges as
   

   ln  ⎛ ⎝ n2 − 1 n2 + 1 ⎞ ⎠ = ln ⎛ ⎝ 1 − 1 n2 ⎞ ⎠ − ln ⎛ ⎝ 1 + 1 n2 ⎞ ⎠  ∼  − 2 n2 .

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.

   (25)

   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,

   (26)

   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).

   (27)

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

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

   (28)

   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 ⎞ ⎠ . (29)

      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). (30)

      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. (31)
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 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,

   (32)

   or
   
   

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

   (33)

   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) . (34)

   Equation (33) can be solved for x as
   
   

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

   (35)

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

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

   (36)

   (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 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .

   (37)

   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 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .

   (38)

   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: