Subsections

• Newton-Cotes Integration of Equations
  • The Trapezoidal rule
  • Simpson's rule
  
  
• Intergrations of Equations
  • Romberg integration
  • Gauss Quadrature
  • Improper integrals
  
  
• Numerical Differentiation
  • High-accuracy differentiation formulas
  • Richardson extrapolation
  • Derivatives of unequally spaced data
  • Numerical integration/differentiation formulas with libraties and packages
  
  
• Engineering Applications: Numerical Integration and Differentiation

Numerical Differentiation and Integration

The derivative represents the rate of cchange of a dependent variable with respect to an independent variable.

$$\frac{dy}{dx} = \lim_{\Delta x \rightarrow 0} \frac{f(x_i+\Delta x)-f(x_i)}{\Delta x}$$ (6.1)

The integration means the total value, or summation, of $f(x)dx$ over the range $x=a$ to $b$.

$$I=\int^a_b f(x)dx$$ (6.2)

Newton-Cotes Integration of Equations

Newton-Gregory forward polynomial : If the $x$-values are evenly spaced, instead of using divided difference, ``ordinary differences'' are more useful; the differences in $f$-values are not divided by the differences in $x$-values.

$$\begin{split}P_n(x_s) =& f_0 + s \Delta f_0 + \frac{s(s-1)}{2... ... &+ \frac{s(s-1)\cdots(s-n+1)}{n\!} \Delta^n f_0 \end{split}$$ (6.3)

where $s=(x-x_0)/h$, with $h=\Delta x$, the uniform spacing in $x$-values.

The Newton-Cotes formulas are based on the strategy of replacing a complicated function or tabulated data with an approximating function that is easy to integrate:

$$I=\int^a_b f(x)dx \simeq \int^a_b P_n(x)dx$$ (6.4)

where $P_n(x)$ is the Newton-Gregory interpolating polynomial. For $n=1$

$$\begin{split}\int^a_b f(x) dx \simeq & \int^a_b P_1(x)dx = \i... ...c{1}{2} \Delta f_0 \right) = \frac{h}{2}(f_0 + f_1) \end{split}$$ (6.5)

For $n=2$

$$\begin{split}\int^a_b f(x) dx &\simeq \int^a_b \left( f_0 + s... ...0 \right) ds & = \frac{h}{3}(f_0 + 4 f_1 + f_2) \end{split}$$ (6.6)

See the figure 21.1 in the textbook.

The Trapezoidal rule

The trapezoidal rule is the first of the Newton-Cotes closed integration formulas

$$I=\int^a_b f(x)dx \simeq \int^a_b f_1(x)dx$$ (6.7)

where

$$f_1(x) = f(a) + \frac{f(b)-f(a)}{b-a}(x-a)$$ (6.8)

The result of integration is

$$I = (b-a) \frac{f(b)-f(a)}{2}$$ (6.9)

which is called as trapezoidal rule.

One way to improve the accuracy of the trapezoidal rule is to divide the integration interval from $a$ to $b$ into a number of segments and apply the method to each segment. The width of segments

$$h=\frac{b-a}{n}$$ (6.10)

The integration is

$$I=\frac{h}{2}\left[ f(x_0) + 2 \sum^{n-1}_{i=1} f(x_i) + f(x_n) \right]$$ (6.11)

Simpson's rule

Another way to obtain a more accurate estimate of an integral is to use higher-order polynomial to connect the points.

Simpson's 1/3 rule : use a second-order polynomial

$$I=\int^b_a f(x)dx \simeq \int^b_a f_2(x)dx$$ (6.12)

Simpson's 1/3 rule is

$$I \simeq \frac{h}{3} \left[ f(x_0) + 4f(x_1) + f(x_2) \right]$$ (6.13)

The label ``1/3'' stems from the fact that $h$ is divided by 3.

Simpson's 3/8 rule : use a third-order Lagrange polynomial

$$I=\int^b_a f(x)dx \simeq \int^b_a f_3(x)dx$$ (6.14)

Simpson's 3/8 rule is

$$I \simeq \frac{3h}{8} \left[ f(x_0) + 3f(x_1) + 3f(x_2) + f(x_3) \right]$$ (6.15)

See the figure 21.11 in the textbook.

Intergrations of Equations

Romberg integration

Richardson's extrapolation : use two estimates of an integral to compute a third. It improves the results of numerical integration on the basis of the integral estimate themselves.

Two separate estimate using step sizes of $h_1$ and $h_2$

$$I(h_1) + E(h_1) = I(h_2) + E(h_1)$$ (6.16)

The error of the multiple-application trapezoidal rule is

$$E \simeq -\frac{b-a}{12}h^2 \bar{f}''$$ (6.17)

Assume that $\bar{f}''$ is constant regardless of step size

$$\frac{E(h_1)}{E(h_2)} \simeq \frac{h^2_1}{h^2_2}$$ (6.18)

Rearranage the above equation

$$E(h_1) \simeq E(h_2) \left(\frac{h_1}{h_2}\right)^2$$ (6.19)

which can be substituted into eq. (6.16)

$$I(h_1) + E(h_2) \left(\frac{h_1}{h_2}\right)^2 \simeq I(h_2) + E(h_1)$$ (6.20)

which can be solved for

$$E(h_2) \simeq \frac{I(h_1)-I(h_2)}{1-(h_1/h_2)^2}$$ (6.21)

Thus, we have developed an estimate of the truncation error in terms of the integral estimates and their step sizes. This estimate can then be substituted into

$$I = I(h_2)+E(h_2)$$ (6.22)

to yield an improved estimate of the integral:

$$I \simeq I(h_2) + \frac{1}{(h_1/h_2)^2 - 1}\left[ I(h_2) - I(h_1) \right]$$ (6.23)

For the special case where the interval is halved $(h_2 = h_1 /2)$

$$I \simeq I(h_2) + \frac{1}{2^2 - 1}\left[ I(h_2) - I(h_1) \right]$$ (6.24)

or

$$I \simeq \frac{4}{3} I(h_2) - \frac{1}{3}I(h_1)$$ (6.25)

The Romberg integration algorithm

$$I_{j,k} \simeq \frac{4^{k-1} I_{j+1,k-1} - I_{j,k-1}}{4^{k-1} - 1}$$ (6.26)

where $I_{j+1,k+1}$ and $I_{j,k-1}$ are the more and less accurate integral and $I{j,k}$ is the improved integral.

Gauss Quadrature

• Trapezoidal rule : two parameters model

  $$I \simeq c_0 f(a) + c_1 f(b)$$ (6.27)

  
  
  where the $c$'s are the unknown parameters.
• Gauss Quadrature : four parameters model

  $$I \simeq c_0 f(a) + c_1 f(b)$$ (6.28)

  
  
  where the $c$'s, $f(a)$, $f(b)$ are the unknown parameters.

The trapezoidal rule's formula can be derived from another point of view, the method of undetermined coefficients. Because the trapezoidal rule is a two parameters model, we need two relationships that connect two parameters.

$$c_0 + c_1 = \int^{(b-a)/2}_{-(b-a)/2} 1 dx$$ (6.29)

and

$$-c_0 \frac{b-a}{2} + c_1 \frac{b-a}{2} = \int^{(b-a)/2}_{-(b-a)/2} x dx$$ (6.30)

$$I \simeq (b-a) \frac{f(a)+f(b)}{2}$$ (6.31)

The trapezoidal rule must pass through the end point and results in a large error. But suppose that the constraints of fixed base points was removed and we were freely evaluate the area under a straight line joining any two points on the curve. See the figure 22.5 to figure out the differences.

The object of Gauss quadrature is to determine the coefficients of an equation of the form

$$I \simeq c_0 f(a) + c_1 f(b)$$ (6.32)

with assuming that eq. (6.32) fit the integral of a constant, a linear, a parabolic, and a cubic function

$$c_0 f(x_0) + c_1 f(x_1) = \int^1_{-1} dx = 2$$ (6.33)

$$c_0 f(x_0) + c_1 f(x_1) = \int^1_{-1} x dx = 0$$ (6.34)

$$c_0 f(x_0) + c_1 f(x_1) = \int^1_{-1} x^2 dx = \frac{2}{3}$$ (6.35)

$$c_0 f(x_0) + c_1 f(x_1) = \int^1_{-1} x^3 dx = 0$$ (6.36)

These relationships yield the two-point Gauss-Legendre formula

$$I \simeq f(-\frac{1}{\sqrt{3}}) + f(\frac{1}{\sqrt{3}})$$ (6.37)

Because Gauss quadrature requires function evalutions at nonuniformly spaced points within the integration interval, it is not appropriate for cases where the function is unknown.

Improper integrals

Improper intergral that is one with a lower limit of $-\infty$ or an upper limit of $+\infty$, usually can be evaluated by makeing a change of variable that transforms the infinite range to one that is finite.

Numerical Differentiation

Improve derivative estimates

• Decrease the step size
• Use a higher-order formula
• Combine two derivative estimates to compute more accurate approximation

High-accuracy differentiation formulas

The forward Taylor series expansion can be written as

$$f(x_{i+1}) = f(x_i) + f'(x_i)h + \frac{f''(x_i)}{2}h^2 + \cdots$$ (6.38)

which can be solved for

$$f'(x_i) = \frac{f(x_{i+1})-f(x_i)}{h} - \frac{f''(x_i)}{2}h^2 + O(h^2)$$ (6.39)

If we truncate the second- and higher-derivative terms

$$f'(x_i) = \frac{f(x_{i+1}-f(x_i)}{h}$$ (6.40)

The accuracy of the above equation depend on the step size $h$.

In contrast to this approach, substitue the second-derivative term

$$f''(x_i) = \frac{f(x_{i+2}) - 2f(x_{i+1}) + f(x_i)}{h^2} + O(h)$$ (6.41)

into eq. (6.39) to yield

$$f'(x_i) = \frac{f(x_{i+1})-f(x_i)}{h} - \frac{f(x_{i+2}) - 2f(x_{i+1}) + f(x_i)}{2h^2}h+ O(h^2)$$ (6.42)

or, by collecting terms,

$$f'(x_i) = \frac{-f(x_{i+2}) + 4f(x_{i+1}) - 3f(x_i)}{2h} + O(h^2)$$ (6.43)

Notice that inclusion of the second-derivative term has improved the accuracy to $O(h^2)$.

Richardson extrapolation

In similar with Richardson extrapolation for integral, an estimate for derivatives can be written as

$$D \simeq \frac{4}{3} D(h_2) - \frac{1}{3} D(h_1)$$ (6.44)

For centered difference approximations with $O(h^2)$, the application of this formula will yield a new derivative estimate of $O(h^4)$.

Derivatives of unequally spaced data

One way to handle nonequispaceddata is to fit a second-order Lagrange interpolating polynomial to each set of three adjacent points.

Numerical integration/differentiation formulas with libraties and packages

Various subroutines and functions are exist to solve integral and derivative problems in Matlab and IMSL.

Engineering Applications: Numerical Integration and Differentiation