Skip to content
PostNumerical Computation / Lecture

05-NC-Numerical Differentiation and Integration

2025-12-26
Back to Blog

Numerical Differentiation and Integration

Numerical Differentiation

The derivative represents the rate of change of a dependent variable y or f(x) with respect to an independent variable x. (derivative = f(x))


There are three basic ways to estimate derivatives: Forward, Backward and Center finite divided difference All based on the Taylor Series

f(x+h)=f(x)+f(x)h+f(x)2!h2+

Another form of this result is the following, where c is some point between x and x+h.

f(x+h)=f(x)+f(x)h+f(c)2!h2

Rearranging this last equation gives

f(x)=f(x+h)f(x)hf(c)2h

We will estimate the derivative by dropping the last term. If we can bound the second derivative by some constant, then this last term is O(h) and we get

f(x)=f(x+h)f(x)h+O(h)

Thus,

f(x)f(x+h)f(x)h

This formula is called the first forward divide difference formula and the error in this formula is of order O(h).

The error in this approximating is due to truncation of the last term (the second derivative term).

The truncation error can be made small by making h small. However, as h gets smaller, precision will be lost in this equation due to subtractive cancellation.

The error in calculation for small h is called round-off error. Generally the round-off error will increase as h decreases.

EXAMPLE:

use first forward divide difference formula with h=0.1 to approximate the derivative of f(x)=1/x at x=2.

SOLUTION:

f(x)f(x+h)f(x)h=12.1120.10.2381.

the error is

0.2381(122)=0.0119.

Compare this to the error predicted by the formula, which is f(c)h/2 for some c between 2 and 2.1. Since f(x)=2x3, the error must be between (0.1)230.0125 and (0.1)(2.1)30.0108, which is consistent with our result. However, this information is usually not available.


More precise solution

Consider the Taylor Series

f(xh)=f(x)f(x)h+f(c)2!h2

Rearranging this last equation gives

f(x)=f(x)f(xh)h+f(c)2h

Or,

f(x)=f(x)f(xh)h+O(h)

This is called the backward divide difference formula for f(x).

Consider the Taylor series expansions for the forward and backward approximations, extended to the degree 3 terms:

f(x+h)=f(x)+f(x)h+f(x)2!h2+f(c1)3!h3f(xh)=f(x)f(x)h+f(x)2!h2f(c2)3!h3

If we subtract these two and solve for f(x) we get:

f(x)=f(x+h)f(xh)2h[f(c1)+f(c2)2]h26

THEOREM Generalized intermediate value theorem

Let f be a continuous function on the interval [a,b]. Let x1,,xn be points in [a,b], and a1,,an>0. Then there exists a number c between a and b such that

(a1++an)f(c)=a1f(x1)++anf(xn).

Based on Generalized intermediate value theorem, we have

f(x)=f(x+h)f(xh)2hh26f(c),

where xh<c<x+h.

f(x)=f(x+h)f(xh)2h+O(h2)

This is called the centered divide difference formula for f(x).

EXAMPLE

use centered divide difference formula with h=0.1 to approximate the derivative of f(x)=1/x at x=2.

SOLUTION:

f(x)f(x+h)f(xh)2h=12.111.90.20.2506.

the error is

0.2506(122)=0.0006.

Second derivative

Consider the Taylor series expansions for the forward and backward approximations, extended to the degree 4 terms:

f(x+h)=f(x)+f(x)h+f(x)2!h2+f(x)3!h3+f(4)(x)4!h4+f(xh)=f(x)f(x)h+f(x)2!h2f(x)3!h3+f(4)(x)4!h4

If we add these two and solve for f(x) we get:

f(x)=f(x+h)2f(x)+f(xh)h2+O(h2)

Try to Improve

The forward difference formula for f(x) had error O(h), while the centered difference formula had error O(h2). Can we do better?

Consider again the Taylor expansions:

f(x+h)=f(x)+f(x)h+f(x)2!h2+f(x)3!h3+f(4)(x)4!h4+f(5)(x)5!h5+f(xh)=f(x)f(x)h+f(x)2!h2f(x)3!h3+f(4)(x)4!h4f(5)(x)5!h5+

Let ϕ(h)=f(x+h)f(xh)2h Then,

ϕ(h)=f(x)+a2h2+a4h4+a6h6+

where a2=f(x)3, a4=f(5)(x)60, ... which can just be considered as some constant coefficients.

Also, ϕ(h2)=f(x)+14a2h2+116a4h4+164a6h6+

We can use ϕ(h) and ϕ(h2) to get:

ϕ(h)4ϕ(h2)=3f(x)+34a4h4+1516a6h6+6364a8h8+4ϕ(h2)ϕ(h)3=f(x)14a4h4516a6h62164a8h8+

Note: This approximation for f(x) has error of O(h4)!!

f(x)4ϕ(h2)ϕ(h)3f(x)4ϕ(h2)ϕ(h)3=[4f(x+h2)f(xh2)hf(x+h)f(xh)2h]3=f(xh)8f(xh2)+8f(x+h2)f(x+h)6h

Richardson Extrapolation

We can use the Richardson method to increase the accuracy of numerical estimates to any series-based quantity. Suppose we want to calculate some quantity L and have found, through theory, an approximation to L:

ϕ(h)=L+k=1a2kh2k

Let D(n,0)=ϕ(h2n), define:

D(n,m)=4mD(n,m1)D(n1,m1)4m1D(n,0)=ϕ(h2n)

THEOREM: (Richardson Extrapolation).

There are constants ak,m such that

D(n,m)=L+k=m+1ak,m(h2n)2k(0mn).

COROLLARY:

D(n,n)=L+O(h2(n+1))

D(n,n) Can be found through a pyramid table:

D(0,0)D(1,0)D(1,1)D(2,0)D(2,1)D(2,2)D(n,0)D(n,1)D(n,2)D(n,n)

EXAMPLE:

Consider f(x)=arctan(x). Suppose we want to find f(2). Let ϕ(h)=12h[f(2+h)f(2h)], and start with h=0.01.

Then, we compute D(n,m) in a pyramid fashion (as we did for Newton’s divided differences). The first column is just D(n,0)=ϕ(h2n)

D(n,m)012
00.333339506181068
10.3333348765437230.333333333331274
20.3333337191358200.3333333333331860.333333333333313

Numerical Integration

Trapezoidal rule

The need often arises for evaluating the definite integral of a function that has no explicit antiderivative or whose antiderivative is not easy to obtain. In this section, we will discuss some numerical method to approximating abf(x)dx.

Let f(x) be a function with a continuous second derivative, defined one the interval [x0,x1], as shown in the figure. Denote the corresponding function values by y0=f(x0) and y1=f(x1). Consider the degree 1 interpolating polynomial P1(x) through (x0,y0) and (x1,y1). Using the Lagrange formulation, we find that the interpolating polynomial with error term is

f(x)=y0xx1x0x1+y1xx0x1x0+(xx0)(xx1)2!f(cx)=P(x)+E(x).

cx depends continuously on x.

Integrating both sides on the interval of interest [x0,x1] yields

x0x1f(x)dx=x0x1P(x)dx+x0x1E(x)dx.

Computing the first integral gives

x0x1P(x)dx=y0x0x1xx1x0x1dx+y1x0x1xx0x1x0dx=y0h2+y1h2=hy0+y12,

where we have defined h=x1x0 to be the interval length.

The error term is

x0x1E(x)dx=12!x0x1(xx0)(xx1)f(c(x))dx=f(c)2x0x1(xx0)(xx1)dx=f(c)20hu(uh)du=h312f(c).x0x1f(x)dx=h2(y0+y1)h312f(c)

where h=x1x0 and c is between x0 and x1.

Simpson Rule

We replace degree 1 interpolant by parabola.

f(x)=y0(xx1)(xx2)(x0x1)(x0x2)+y1(xx0)(xx2)(x1x0)(x1x2)+y2(xx0)(xx1)(x2x0)(x2x1)+(xx0)(xx1)(xx2)3!f(cx)=P(x)+E(x).

Integrating gives

x0x2f(x)dx=x0x2P(x)dx+x0x2E(x)dx,

where

x0x2P(x)dx=y0x0x2(xx1)(xx2)(x0x1)(x0x2)dx+y1x0x2(xx0)(xx2)(x1x0)(x1x2)dx+y2x0x2(xx0)(xx1)(x2x0)(x2x1)dx=y0h3+y14h3+y2h3.

where h=x2x1=x1x0.

The error term can be computed (proof omitted) as

x0x2E(x)dx=h590f(4)(c)

for some c in the interval [x0,x2], provided that f(4)(c) exists and is continuous. Concluding the derivation yields Simpson’s rule

Simpson’s Rule

x0x2f(x)dx=h3(y0+4y1+y2)h590f(4)(c),

where h=x1x0=x2x1 and c is between x0 and x2.


EXAMPLE

Apply the Trapezoid rule and Simpson’s rule to approximate

12lnxdx,

and find an upper bound for the error in your approximations.

SOLUTION

12lnxdxh2(y0+y1)=12(ln1+ln2)=ln220.3466.

The error for the Trapezoid rule is h312f(c), where 1<c<2. Since f(x)=1/x2, the magnitude of the error is at most

1312c21120.0834

In other words,

12lnxdx=0.3466±0.083412lnxdxh3(y0+4y1+y2)=0.53(ln1+4ln32+ln2)0.3858.

The error for the Simpson’s rule is h590f(4)(c), where 1<c<2. Since f(4)(x)=6/x4, the magnitude of the error is at most

6(0.5)590c46(0.5)590=14800.0021.

In other words,

12lnxdx=0.3858±0.0021.

The integral can be computed exactly by using integration by parts:

12lnxdx=xlnx|12121dx=2ln21ln110.386294.

The Trapezoid rule and Simpson’s rule error bounds are both consistent with this result.

Composite Trapezoidal Rule

Consider an evenly spaced grid

a=x0<x1<x2<<xm2<xm1<xm=b

Along the horizontal axis, where h=xi+1xi for each i as shown in figure. On each subinterval, we make the approximation with error term

xixi+1f(x)dx=h2(f(xi)+f(xi+1))h312f(ci),

Assuming that f is continuous, adding up over all subintervals yields

abf(x)dx=h2(f(a)+f(b)+2i=1m1f(xi))i=0m1h312f(ci),

The error term can be written

h312i=0m1f(ci)=h312mf(c),

Since mh=(ba), the error term is (ba)h212f(c)=O(h2).

Composite Trapezoid Rule

abf(x)dx=h2(y0+ym+2i=1m1yi)(ba)h212f(c),

where h=(ba)/m and c is between a and b.

Consider an evenly spaced grid

a=x0<x1<x2<<x2m2<x2m1<x2m=b

Along the horizontal axis, where h=xi+1xi for each i as shown in figure. On each subinterval by the interpolating parabola fit at x2i,x2i+1, and x2i+2, we make the approximation with error term

x2ix2i+2f(x)dx=h3(f(x2i)+4f(x2i+1)+f(x2i+2))h590f(4)(ci).

Composite Simpson’s Rule

Consider an evenly spaced grid

a=x0<x1<x2<<x2m2<x2m1<x2m=b

Along the horizontal axis, where h=xi+1xi for each i as shown in figure. On each subinterval by the interpolating parabola fit at x2i,x2i+1, and x2i+2, we make the approximation with error term

x2ix2i+2f(x)dx=h3(f(x2i)+4f(x2i+1)+f(x2i+2))h590f(4)(ci).

This time, the overlapping is over even-numbered xj only. Adding up over all subintervals yields

abf(x)dx=h3(f(a)+f(b)+4i=1mf(x2i1)+2i=1m1f(x2i))i=0m1h590f(4)(ci).

The error term can be written

h590i=0m1f(4)(ci)=h590mf(4)(c),

Since m2h=(ba), the error term is (ba)h4180f(4)(c)=O(h4).

Composite Simpson’s Rule

abf(x)dx=h3(y0+y2m+4i=1my2i1+2i=1m1y2i)(ba)h4180f(4)(c),

where c is between a and b.

EXAMPLE

Carry out four-panel approximations of

12lnxdx,

Using the composite Trapezoid rule and composite Simpson’s rule.

SOLUTION For composite Trapezoid rule on [1,2], four panels means that h=1/4.

12lnxdx18[y0+y4+2i=13yi]=18[ln1+ln2+2(ln54+ln64+ln74)]0.3837

The error is at most

(ba)h212|f(c)|=1161121c21(16)(12)(12)=11920.0052

For four-panel Simpson’s rule set h=1/8.

12lnxdx1813[y0+y8+4i=14y2i1+2i=13y2i]=124[ln1+ln2+4(ln98+ln118+ln138+ln158)+2(ln54+ln64+ln74)]0.386292

The error is at most

(ba)h4180|f(4)(c)|=(18)411806c41(84)(180)(14)0.000008

EXAMPLE

Find the number of panels m necessary for the composite Simpson’s rule to approximate

0πsin2xdx,

within six correct decimal places.

SOLUTION

We require the error to satisfy

(π0)h4180|f(4)(c)|<0.5×106.

Since the fourth derivative of sin2x is 8cos2x, we need

πh41808<0.5×106,

or h<0.0435. therefore, m=ceil(π2h)=37 panels will be sufficient.

Romberg integration

In this section, let’s use the Richardson Extrapolation technique to improve integral approximations.

Consider the Trapezoidal Rule using 2n subintervals of [a,b]. Then, h=ba2n and the approximation will be

ϕ(n)=12ba2ni=02n1[f(xi)+f(xi+1)]

We know that this approximation has error O(h2). Thus,

ϕ(n)=abf(x)dx+a2h2+a4h4+a6h6+

This is exactly the kind of expression we had for Richardson Extrapolation. Let

R(n,0)=ϕ(n)R(n,m)=4mR(n,m1)R(n1,m1)4m1

Compute the table of values:

R(0,0)R(1,0)R(1,1)R(2,0)R(2,1)R(2,2)R(n,0)R(n,1)R(n,2)R(n,n)

This is called Romberg Integration

EXAMPLE

Approximate 12lnxdx using Romberg integration up to n=3.

SOLUEION

R(0,0)=ϕ(0)=12×(21)×(ln1+ln2)=0.3466R(1,0)=ϕ(1)=12×212×(ln1+2ln32+ln2)=0.3760R(2,0)=ϕ(2)=12×2122×(ln1+2ln54+2ln32+2ln74+ln2)=0.3837R(1,1)=4R(1,0)R(0,0)41=0.3858R(2,1)=4R(2,0)R(1,0)41=0.3863R(2,2)=42R(2,1)R(1,1)421=0.3863
TrapezoidSimpson's
h=10.3466
h=0.500.37600.3858
h=0.250.38370.38630.3863
error3.97×1034.94×1046×106
O(h2)O(h4)O(h6)