Logo
Overview
Numerical Analysis 2

Numerical Analysis 2

Aimless Equation Aimless Equation
May 17, 2026
12 min read

Today I feel bored, so let’s go over some mathematics. We’ll review weeks 4−64-6 even though it’s way past week 4 because matrices are really tedious to typeset in LaTeX\LaTeX so I don’t know if I even want to write about them yet to be honest.

Approximating definite integrals

There’s probably (is) an infinite number of ways to approximate definite integrals. We’ll go over a few of them.

Gaussian quadrature

First, we want to find some way to approximate

∫abf(x) dx.\begin{equation} \int_a^b f(x)\,dx. \end{equation}

We can do this by considering the integration rule

∫abf(x) dx≈A0f(x0)+A1f(x1)+⋯+Anf(xn)=∑i=0nAif(xi),\begin{equation} \int_a^b f(x)\,dx \approx A_0f(x_0) + A_1f(x_1) + \cdots + A_nf(x_n) =\sum_{i=0}^n A_if(x_i), \end{equation}

where x0,x1,…,xnx_0, x_1, \dots, x_n are the nodes in the interval [a,b][a,b] and A0,A1,…,AnA_0, A_1, \dots, A_n are the coefficients. Since there are 2n+22n+2 unknowns, we can make eqref(2) exact for polynomials of degree at most 2n+12n+1 i.e., f(x)=xkf(x) = x^k for k=0,1,…,2n+1k = 0, 1, \dots, 2n+1. This is known as Gaussian quadrature. Now let us look at a few different Gaussian quadrature rules.

Deriving the two-point Gaussian rule

We’ll consider the interval [−d,d][-d,d] and let

∫−ddf(x) dx≈A0f(x0)+A1f(x1).\begin{equation} \int_{-d}^d f(x)\,dx \approx A_0f(x_0) + A_1f(x_1). \end{equation}

eqref(3) gives us the exact result whenever f(x)f(x) is a polynomial of degree at most 3. That is, for f(x)=a0+a1x+a2x2+a3x3f(x)=a_0 + a_1x + a_2x^2 + a_3x^3, for some constants a0,a1,a2,a3a_0, a_1, a_2, a_3. Why you ask? Integrating said constants will tell us.

∫(a0+a1x+a2x2+a3x3) dx=a0∫1 dx+a1∫x dx+a2∫x2 dx+a3∫x3 dx.\begin{equation} \begin{aligned} \int (a_0 + a_1x + a_2x^2 + a_3x^3)\,dx ={}& a_0\int 1\,dx \\ &+ a_1\int x\,dx \\ &+ a_2\int x^2\,dx \\ &+ a_3\int x^3\,dx. \end{aligned} \end{equation}

So eqref(3) will give us the exact results for f(x)=1,x,x2,x3f(x)=1,x,x^2,x^3. It follows that

f(x)=1⟹A0+A1=∫−dd1 dx=2d.\begin{equation} f(x)=1 \Longrightarrow A_0 + A_1 = \int_{-d}^{d} 1\,dx = 2d. \end{equation} f(x)=x⟹A0x0+A1x1=∫−ddx dx=0.\begin{equation} f(x)=x \Longrightarrow A_0x_0 + A_1x_1 = \int_{-d}^{d} x\,dx = 0. \end{equation} f(x)=x2⟹A0x02+A1x12=∫−ddx2 dx=23d3.\begin{equation} f(x)=x^2 \Longrightarrow A_0x_0^2 + A_1x_1^2 = \int_{-d}^{d} x^2\,dx = \frac{2}{3}d^3. \end{equation} f(x)=x3⟹A0x03+A1x13=∫−ddx3 dx=0.\begin{equation} f(x)=x^3 \Longrightarrow A_0x_0^3 + A_1x_1^3 = \int_{-d}^{d} x^3\,dx = 0. \end{equation}

After solving this system of nonlinear equations, notice that A0=A1=dA_0 = A_1 = d and the nodes are symmetric about 00 with magnitude d3\frac{d}{\sqrt{3}}. Without loss of generality, we can take x0=−d3x_0 = -\frac{d}{\sqrt{3}} and x1=d3x_1 = \frac{d}{\sqrt{3}} and the formula becomes

∫−ddf(x) dx≈d[f(−d3)+f(d3)].\begin{equation} \int_{-d}^d f(x)\,dx \approx d\left[ f\left(-\frac{d}{\sqrt{3}}\right) + f\left(\frac{d}{\sqrt{3}}\right) \right]. \end{equation}

Two-point Gaussian rule on [−1,1][-1,1]

Generally speaking, we use the interval [−1,1][-1,1] instead of [−d,d][-d,d], so eqref(9) becomes

∫−11f(x) dx≈f(−13)+f(13).\begin{equation} \int_{-1}^1 f(x)\,dx \approx f\left(-\frac{1}{\sqrt{3}}\right) + f\left(\frac{1}{\sqrt{3}}\right). \end{equation}

Mapping the rule to an arbitrary interval [a,b][a,b]

If you decide to use an arbitrary interval [a,b][a,b], we can do this fairly easily by mapping x∈[−1,1]x\in[-1,1] to t∈[a,b]t\in[a,b] from the transformation function

t=12[(b−a)x+(a+b)].\begin{equation} t = \frac{1}{2}\left[(b-a)x + (a+b)\right]. \end{equation}

eqref(11) forces t=at=a for x=−1x=-1 and t=bt=b for x=1x=1, and dtdx=12(b−a)\frac{dt}{dx} = \frac{1}{2}(b-a). Thus, we have

∫abf(t) dt=(b−a2)∫−11f(12[(b−a)x+(a+b)]) dx.\int_a^b f(t)\,dt = \left(\frac{b-a}{2}\right)\int_{-1}^1 f\left(\frac{1}{2}\left[(b-a)x + (a+b)\right]\right)\,dx.

Since

x0=−13,x1=13,x_0=-\frac{1}{\sqrt{3}}, \qquad x_1=\frac{1}{\sqrt{3}},

the points in [a,b][a,b] are

τ0=a+b2−b−a23,τ1=a+b2+b−a23.\tau_0=\frac{a+b}{2}-\frac{b-a}{2\sqrt{3}}, \qquad \tau_1=\frac{a+b}{2}+\frac{b-a}{2\sqrt{3}}.

So, on the interval [a,b][a,b], the rule becomes

∫abf(t) dt≈(b−a2)[f(a+b2−b−a23)+f(a+b2+b−a23)].\int_a^b f(t)\,dt \approx \left(\frac{b-a}{2}\right)\left[ f\left(\frac{a+b}{2}-\frac{b-a}{2\sqrt{3}}\right) + f\left(\frac{a+b}{2}+\frac{b-a}{2\sqrt{3}}\right) \right].

Let’s go through a quick example.

Example (Approximating an integral with two-point Gaussian quadrature)

Say we want to approximate the following integral and find the absolute error:

∫01ex dx.\int_0^1 e^x\,dx.
Solution

Using the two-point Gaussian rule on [a,b][a,b],

Q2=(b−a2)[f(a+b2−b−a23)+f(a+b2+b−a23)].Q_2=\left(\frac{b-a}{2}\right)\left[ f\left(\frac{a+b}{2}-\frac{b-a}{2\sqrt{3}}\right) + f\left(\frac{a+b}{2}+\frac{b-a}{2\sqrt{3}}\right) \right].

Here f(t)=etf(t)=e^t, a=0a=0, and b=1b=1. Hence,

b−a2=12,a+b2=12,b−a23=123.\frac{b-a}{2}=\frac{1}{2}, \qquad \frac{a+b}{2}=\frac{1}{2}, \qquad \frac{b-a}{2\sqrt{3}}=\frac{1}{2\sqrt{3}}.

So the two nodes are

τ0=12−123,τ1=12+123.\tau_0=\frac{1}{2}-\frac{1}{2\sqrt{3}}, \qquad \tau_1=\frac{1}{2}+\frac{1}{2\sqrt{3}}.

Plugging these into the rule gives

Q2=12[e12−123+e12+123]≈1.717896378.Q_2=\frac{1}{2}\left[e^{\frac{1}{2}-\frac{1}{2\sqrt{3}}} + e^{\frac{1}{2}+\frac{1}{2\sqrt{3}}}\right]\approx 1.717896378.

The exact value is

∫01ex dx=e−1.\int_0^1 e^x\,dx = e-1.

Since e−1≈1.718281828e-1 \approx 1.718281828, the absolute error is

∣Error∣=∣(e−1)−Q2∣≈∣1.718281828−1.717896378∣=3.85450×10−4.|\text{Error}| = \left|(e-1)-Q_2\right| \approx \left|1.718281828-1.717896378\right| =3.85450\times 10^{-4}.

General Gaussian error term

Proposition (General Gaussian error term)

For the error formula below, assume ff has a continuous (2n+2)(2n+2)th derivative on [a,b][a,b]. Then the error from approximating ∫abf(x) dx\int_a^b f(x)\,dx with the (n+1)(n+1)-point Gaussian rule is

E=∫abR2n(x) dx=∫abf(x) dx−[A0f(x0)+A1f(x1)+⋯+Anf(xn)]=∫abf(x) dx−∑i=0nAif(xi)=f(2n+2)(ξ)(2n+2)!∫ab(x−x0)2(x−x1)2⋯(x−xn)2 dx=f(2n+2)(ξ)(2n+2)!∫ab∏i=0n(x−xi)2 dx=f(2n+2)(ξ)(2n+2)!∫abωn2(x) dx,\begin{equation} \begin{aligned} E &= \int_a^b R_{2n}(x)\,dx \\ &= \int_a^b f(x)\,dx - \left[A_0f(x_0)+A_1f(x_1)+\cdots+A_nf(x_n)\right] \\ &= \int_a^b f(x)\,dx - \sum_{i=0}^n A_if(x_i) \\ &= \frac{f^{(2n+2)}(\xi)}{(2n+2)!}\int_a^b (x-x_0)^2(x-x_1)^2\cdots(x-x_n)^2\,dx \\ &= \frac{f^{(2n+2)}(\xi)}{(2n+2)!}\int_a^b \prod_{i=0}^n (x-x_i)^2\,dx \\ &= \frac{f^{(2n+2)}(\xi)}{(2n+2)!}\int_a^b \omega_n^2(x)\,dx, \end{aligned} \end{equation}

for some ξ∈(a,b)\xi\in(a,b). Here R2n(x)R_{2n}(x) is the interpolation remainder, and

ωn(x)=(x−x0)(x−x1)⋯(x−xn)=∏i=0n(x−xi)\omega_n(x)=(x-x_0)(x-x_1)\cdots(x-x_n) =\prod_{i=0}^n (x-x_i)

is the node polynomial. Let

M=max⁡x∈[a,b]∣f(2n+2)(x)∣.M=\max_{x\in[a,b]}\left|f^{(2n+2)}(x)\right|.

Then

∣E∣≤M(2n+2)!∫abωn2(x) dx.\begin{equation} |E|\leq \frac{M}{(2n+2)!}\int_a^b \omega_n^2(x)\,dx. \end{equation}

Two-point error term

Proposition (Two-point Gaussian error term)

For the two-point rule, n=1n=1, so the nodes on [a,b][a,b] are

x0=a+b2−b−a23,x1=a+b2+b−a23.x_0=\frac{a+b}{2}-\frac{b-a}{2\sqrt{3}}, \qquad x_1=\frac{a+b}{2}+\frac{b-a}{2\sqrt{3}}.

Using Proposition 1.3 with n=1n=1, the error becomes

E=f(4)(ξ)4!∫ab(x−a+b2+b−a23)2(x−a+b2−b−a23)2 dx=f(4)(ξ)4!⋅(b−a)5180=(b−a)54320f(4)(ξ),\begin{equation} \begin{aligned} E &= \frac{f^{(4)}(\xi)}{4!} \int_a^b \left(x-\frac{a+b}{2}+\frac{b-a}{2\sqrt{3}}\right)^2 \left(x-\frac{a+b}{2}-\frac{b-a}{2\sqrt{3}}\right)^2 \,dx \\ &= \frac{f^{(4)}(\xi)}{4!}\cdot\frac{(b-a)^5}{180} \\ &= \frac{(b-a)^5}{4320}f^{(4)}(\xi), \end{aligned} \end{equation}

for some ξ∈(a,b)\xi\in(a,b).

Newton-Cotes

Unlike Gaussian quadrature, Newton-Cotes formulas use equally spaced nodes (this will be important later). For a closed Newton-Cotes rule, the endpoints are included, while as you’d expect, an open Newton-Cotes rule excludes the endpoints.

xi=a+ih,h=b−an,i=0,1,…,n,n≥1.\begin{equation} x_i=a+ih, \qquad h=\frac{b-a}{n}, \qquad i=0,1,\dots,n, \qquad n\geq 1. \end{equation}

The general form is still

∫abf(x) dx≈∑i=0nAif(xi),\begin{equation} \int_a^b f(x)\,dx \approx \sum_{i=0}^n A_if(x_i), \end{equation}

We’ll be looking at the Trapezoidal rule and Simpson’s rule which comes from the closed Newton-Cotes formula with n=1n=1 and n=2n=2 respectively. But first, let’s get the general error terms for closed Newton-Cotes formulas out of the way.

Proposition (Closed Newton-Cotes error term)

Suppose

∫abf(x) dx≈∑i=0nAif(xi)\begin{equation} \int_a^b f(x)\,dx \approx \sum_{i=0}^n A_if(x_i) \end{equation}

is the (n+1)(n+1)-point closed Newton-Cotes formula with x0=ax_0=a, xn=bx_n=b, h=b−anh=\frac{b-a}{n}, and n≥1n\geq 1. Define the error by

En=∫abf(x) dx−∑i=0nAif(xi).\begin{equation} E_n=\int_a^b f(x)\,dx-\sum_{i=0}^n A_if(x_i). \end{equation}

Ck([a,b])C^k([a,b]) denotes the set of functions with kk continuous derivatives on [a,b][a,b].

If nn is even and f∈Cn+2([a,b])f\in C^{n+2}([a,b]), then there exists ξ∈(a,b)\xi\in(a,b) such that

En=hn+3f(n+2)(ξ)(n+2)!∫0nt2(t−1)⋯(t−n) dt=hn+3f(n+2)(ξ)(n+2)!∫0nt2∏i=1n(t−i) dt\begin{equation} E_n= \frac{h^{n+3}f^{(n+2)}(\xi)}{(n+2)!} \int_0^n t^2(t-1)\cdots(t-n)\,dt = \frac{h^{n+3}f^{(n+2)}(\xi)}{(n+2)!} \int_0^n t^2\prod_{i=1}^n(t-i)\,dt \end{equation}

If nn is odd and f∈Cn+1([a,b])f\in C^{n+1}([a,b]), then there exists ξ∈(a,b)\xi\in(a,b) such that

En=hn+2f(n+1)(ξ)(n+1)!∫0nt(t−1)⋯(t−n) dt=hn+2f(n+1)(ξ)(n+1)!∫0n∏i=0n(t−i) dt\begin{equation} E_n= \frac{h^{n+2}f^{(n+1)}(\xi)}{(n+1)!} \int_0^n t(t-1)\cdots(t-n)\,dt = \frac{h^{n+2}f^{(n+1)}(\xi)}{(n+1)!} \int_0^n \prod_{i=0}^n(t-i)\,dt \end{equation}
Explanation

I split the cases to make the distinction clear. Specifically the parity of nn. The Trapezoidal rule has n=1n=1, so it uses the odd case in eqref(20). Simpson’s rule has n=2n=2, so it uses the even case in eqref(19). Following the degree we can see that the Trapezoidal rule approximates a curve with a straight line, while Simpson’s rule approximates the curve with a parabola. Both Trapezoidal and Simpson’s rule are from the closed Newton-Cotes family. Now we know that any Newton-Cotes formula uses equally spaced nodes, which includes the Trapezoidal and Simpson’s rules. It’s very rigid in that sense, the data we use to approximate the integral must be equally spaced. For data on a uniform grid it’s nice but we aren’t free to choose better points. Gaussian quadrature instead treats the locations of these sample points as part of the problem. Instead of choosing equally spaced nodes to simplify the problem (since they only have to deal with the weights), the goal is to find the “best” nodes and weights. So the points aren’t arbitrary, but placed in locations that make the formula exact for the largest possible range of polynomial behaviour. To avoid too much ambiguity, we don’t use evenly spaced points/nodes or “more” points to get a better approximation, rather “better” positioned points; and these “better” points are placed where they give the most information about the integral. That is, the roots of the polynomials. Now let’s actually look at the approximation methods.

Trapezoidal rule

Proposition (Trapezoidal rule)

The Trapezoidal rule is the closed Newton-Cotes rule with n=1n=1. Therefore,

x0=a,x1=b,h=b−a.\begin{equation} x_0=a, \qquad x_1=b, \qquad h=b-a. \end{equation}

Approximating ff by (a,f(a))(a,f(a)) and (b,f(b))(b,f(b)) gives

∫abf(x) dx≈b−a2[f(a)+f(b)].\begin{equation} \int_a^b f(x)\,dx \approx \frac{b-a}{2}\left[f(a)+f(b)\right]. \end{equation}

The error term is

ET=−(b−a)312f′′(ξ)\begin{equation} E_T=-\frac{(b-a)^3}{12}f''(\xi) \end{equation}

for some ξ∈(a,b)\xi\in(a,b). Thus,

∣ET∣≤(b−a)312max⁡x∈[a,b]∣f′′(x)∣.\begin{equation} |E_T|\leq \frac{(b-a)^3}{12}\max_{x\in[a,b]}\left|f''(x)\right|. \end{equation}
Example (Trapezoidal approximation)

Approximate

∫01ex dx\int_0^1 e^x\,dx

using the trapezoidal rule, and find an upper bound for the error.

Solution

Here a=0a=0 and b=1b=1. Using eqref(22),

T=12[f(0)+f(1)]=12(1+e)≈1.859140914.\begin{equation} T=\frac{1}{2}\left[f(0)+f(1)\right] =\frac{1}{2}(1+e) \approx 1.859140914. \end{equation}

For f(x)=exf(x)=e^x, we have f′′(x)=exf''(x)=e^x. On [0,1][0,1],

max⁡x∈[0,1]∣f′′(x)∣=e.\begin{equation} \max_{x\in[0,1]}\left|f''(x)\right|=e. \end{equation}

Using eqref(24),

∣ET∣≤(1−0)312e=e12≈2.26523×10−1.\begin{equation} |E_T|\leq \frac{(1-0)^3}{12}e=\frac{e}{12}\approx 2.26523\times 10^{-1}. \end{equation}

Since the exact value is e−1e-1, the actual absolute error is

∣(e−1)−T∣≈1.40859×10−1.\begin{equation} \left|(e-1)-T\right|\approx 1.40859\times 10^{-1}. \end{equation}

Simpson’s rule

Proposition (Simpson's rule)

Simpson’s rule is the closed Newton-Cotes rule with n=2n=2. The three equally spaced nodes are

x0=a,x1=a+b2,x2=b.\begin{equation} x_0=a, \qquad x_1=\frac{a+b}{2}, \qquad x_2=b. \end{equation}

Approximating ff gives

∫abf(x) dx≈b−a6[f(a)+4f(a+b2)+f(b)].\begin{equation} \int_a^b f(x)\,dx \approx \frac{b-a}{6}\left[ f(a)+4f\left(\frac{a+b}{2}\right)+f(b) \right]. \end{equation}

The error term is

ES=−(b−a)52880f(4)(ξ)\begin{equation} E_S=-\frac{(b-a)^5}{2880}f^{(4)}(\xi) \end{equation}

for some ξ∈(a,b)\xi\in(a,b). Thus,

∣ES∣≤(b−a)52880max⁡x∈[a,b]∣f(4)(x)∣.\begin{equation} |E_S|\leq \frac{(b-a)^5}{2880}\max_{x\in[a,b]}\left|f^{(4)}(x)\right|. \end{equation}
Example (Simpson approximation)

Approximate

∫01ex dx\int_0^1 e^x\,dx

using Simpson’s rule, and find an upper bound for the error.

Solution

Here a=0a=0, b=1b=1, and a+b2=12\frac{a+b}{2}=\frac{1}{2}. Using eqref(30),

S=16[f(0)+4f(12)+f(1)]=16[1+4e1/2+e]≈1.718861152.\begin{equation} S=\frac{1}{6}\left[f(0)+4f\left(\frac{1}{2}\right)+f(1)\right] =\frac{1}{6}\left[1+4e^{1/2}+e\right] \approx 1.718861152. \end{equation}

For f(x)=exf(x)=e^x, we have f(4)(x)=exf^{(4)}(x)=e^x. On [0,1][0,1],

max⁡x∈[0,1]∣f(4)(x)∣=e.\begin{equation} \max_{x\in[0,1]}\left|f^{(4)}(x)\right|=e. \end{equation}

Using eqref(32),

∣ES∣≤(1−0)52880e=e2880≈9.43848×10−4.\begin{equation} |E_S|\leq \frac{(1-0)^5}{2880}e=\frac{e}{2880}\approx 9.43848\times 10^{-4}. \end{equation}

Since the exact value is e−1e-1, the actual absolute error is

∣(e−1)−S∣≈5.79323×10−4.\begin{equation} \left|(e-1)-S\right|\approx 5.79323\times 10^{-4}. \end{equation}

Open Newton-Cotes formulas

Open Newton-Cotes formulas still use equally spaced nodes, but they do not use the endpoints. Instead, the endpoints sit one step outside the nodes:

x−1=a,xn+1=b,h=b−an+2.\begin{equation} x_{-1}=a, \qquad x_{n+1}=b, \qquad h=\frac{b-a}{n+2}. \end{equation}

So the nodes used by the rule are

xi=a+(i+1)h,i=0,1,…,n.\begin{equation} x_i=a+(i+1)h, \qquad i=0,1,\dots,n. \end{equation}
Proposition (Open Newton-Cotes error term)

Suppose

∫abf(x) dx≈∑i=0nAif(xi)\begin{equation} \int_a^b f(x)\,dx \approx \sum_{i=0}^n A_if(x_i) \end{equation}

is the (n+1)(n+1)-point open Newton-Cotes formula with x−1=ax_{-1}=a, xn+1=bx_{n+1}=b, and h=b−an+2h=\frac{b-a}{n+2}. We define the error by

En=∫abf(x) dx−∑i=0nAif(xi).\begin{equation} E_n=\int_a^b f(x)\,dx-\sum_{i=0}^n A_if(x_i). \end{equation}

Here Ck([a,b])C^k([a,b]) denotes the set of functions with kk continuous derivatives on [a,b][a,b].

If nn is even and f∈Cn+2([a,b])f\in C^{n+2}([a,b]), then there exists ξ∈(a,b)\xi\in(a,b) such that

En=hn+3f(n+2)(ξ)(n+2)!∫−1n+1t2(t−1)⋯(t−n) dt=hn+3f(n+2)(ξ)(n+2)!∫−1n+1t2∏i=1n(t−i) dt.\begin{equation} E_n= \frac{h^{n+3}f^{(n+2)}(\xi)}{(n+2)!} \int_{-1}^{n+1} t^2(t-1)\cdots(t-n)\,dt = \frac{h^{n+3}f^{(n+2)}(\xi)}{(n+2)!} \int_{-1}^{n+1} t^2\prod_{i=1}^n(t-i)\,dt. \end{equation}

If nn is odd and f∈Cn+1([a,b])f\in C^{n+1}([a,b]), then there exists ξ∈(a,b)\xi\in(a,b) such that

En=hn+2f(n+1)(ξ)(n+1)!∫−1n+1t(t−1)⋯(t−n) dt=hn+2f(n+1)(ξ)(n+1)!∫−1n+1∏i=0n(t−i) dt.\begin{equation} E_n= \frac{h^{n+2}f^{(n+1)}(\xi)}{(n+1)!} \int_{-1}^{n+1} t(t-1)\cdots(t-n)\,dt = \frac{h^{n+2}f^{(n+1)}(\xi)}{(n+1)!} \int_{-1}^{n+1} \prod_{i=0}^n(t-i)\,dt. \end{equation}

As you can see from Proposition 1.5, it’s pretty much exactly the same as the closed case, except for the being open part (and some shifted indexing, spacing, etc.).

Composite formulas

Instead of applying one rule across the whole interval, composite formulas split [a,b][a,b] into smaller subintervals and apply the same rule repeatedly.

Composite Trapezoidal rule

Proposition (Composite trapezoidal rule)

Let [a,b][a,b] be divided into nn equal subintervals, where

h=b−an,xi=a+ih,i=0,1,…,n.\begin{equation} h=\frac{b-a}{n}, \qquad x_i=a+ih, \qquad i=0,1,\dots,n. \end{equation}

Then

∫abf(x) dx≈h2[f(x0)+2∑i=1n−1f(xi)+f(xn)].\begin{equation} \int_a^b f(x)\,dx \approx \frac{h}{2}\left[ f(x_0)+2\sum_{i=1}^{n-1}f(x_i)+f(x_n) \right]. \end{equation}

If M=max⁡x∈[a,b]∣f′′(x)∣M=\max_{x\in[a,b]}\left|f''(x)\right|, then the error satisfies

∣ET∣≤Mnh312.\begin{equation} |E_T|\leq \frac{Mnh^3}{12}. \end{equation}
Example (Composite trapezoidal approximation)

Approximate

∫01ex dx\int_0^1 e^x\,dx

using the composite trapezoidal rule with 66 equal subintervals. Then determine an error bound and compare it with the true error.

Solution

Here a=0a=0, b=1b=1, and n=6n=6, so

h=1−06=16,xi=i6.\begin{equation} h=\frac{1-0}{6}=\frac{1}{6}, \qquad x_i=\frac{i}{6}. \end{equation}
iixix_if(xi)=exif(x_i)=e^{x_i}
00001.0000000001.000000000
1116\frac{1}{6}1.1813604131.181360413
2213\frac{1}{3}1.3956124251.395612425
3312\frac{1}{2}1.6487212711.648721271
4423\frac{2}{3}1.9477340411.947734041
5556\frac{5}{6}2.3009758912.300975891
66112.7182818282.718281828

Using eqref(44),

T6=h2[f(x0)+2∑i=15f(xi)+f(x6)]=112[1+2∑i=15ei/6+e]=112[1+2(1.181360413+1.395612425+1.648721271+1.947734041+2.300975891)+2.718281828]≈1.722257492.\begin{equation} \begin{aligned} T_6 &=\frac{h}{2}\left[ f(x_0)+2\sum_{i=1}^{5}f(x_i)+f(x_6) \right] \\ &=\frac{1}{12}\left[ 1+2\sum_{i=1}^{5}e^{i/6}+e \right] \\ &=\frac{1}{12}\left[ 1+2(1.181360413+1.395612425+1.648721271+1.947734041+2.300975891)+2.718281828 \right] \\ &\approx 1.722257492. \end{aligned} \end{equation}

For f(x)=exf(x)=e^x, we have f′′(x)=exf''(x)=e^x. On [0,1][0,1],

M=max⁡x∈[0,1]∣f′′(x)∣=e.\begin{equation} M=\max_{x\in[0,1]}\left|f''(x)\right|=e. \end{equation}

Using eqref(45),

∣ET∣≤Mnh312=e⋅6⋅(16)312=e432≈6.29232×10−3.\begin{equation} |E_T|\leq \frac{Mnh^3}{12} =\frac{e\cdot 6\cdot \left(\frac{1}{6}\right)^3}{12} =\frac{e}{432} \approx 6.29232\times 10^{-3}. \end{equation}

The exact value is e−1e-1, so the true error is

∣(e−1)−T6∣≈∣1.718281828−1.722257492∣=3.97566×10−3.\begin{equation} \left|(e-1)-T_6\right| \approx \left|1.718281828-1.722257492\right| =3.97566\times 10^{-3}. \end{equation}

Composite Simpson’s rule

Proposition (Composite Simpson's rule)

Let [a,b][a,b] be divided into 2n2n equal subintervals, where

h=b−a2n,xi=a+ih,i=0,1,…,2n.\begin{equation} h=\frac{b-a}{2n}, \qquad x_i=a+ih, \qquad i=0,1,\dots,2n. \end{equation}

Then

∫abf(x) dx≈h3[f(x0)+4∑i=1nf(x2i−1)+2∑i=1n−1f(x2i)+f(x2n)].\begin{equation} \int_a^b f(x)\,dx \approx \frac{h}{3}\left[ f(x_0) +4\sum_{i=1}^{n}f(x_{2i-1}) +2\sum_{i=1}^{n-1}f(x_{2i}) +f(x_{2n}) \right]. \end{equation}

If M=max⁡x∈[a,b]∣f(4)(x)∣M=\max_{x\in[a,b]}\left|f^{(4)}(x)\right|, then the error satisfies

∣ES∣≤Mnh590.\begin{equation} |E_S|\leq \frac{Mnh^5}{90}. \end{equation}
Example (Composite Simpson approximation)

Approximate

∫01ex dx\int_0^1 e^x\,dx

using the composite Simpson’s rule with 66 equal subintervals. Then determine an error bound and compare it with the true error.

Solution

Here 2n=62n=6, so n=3n=3. Also a=0a=0, b=1b=1, and

h=1−06=16,xi=i6.\begin{equation} h=\frac{1-0}{6}=\frac{1}{6}, \qquad x_i=\frac{i}{6}. \end{equation}
iixix_if(xi)=exif(x_i)=e^{x_i}Simpson weight
00001.0000000001.00000000011
1116\frac{1}{6}1.1813604131.18136041344
2213\frac{1}{3}1.3956124251.39561242522
3312\frac{1}{2}1.6487212711.64872127144
4423\frac{2}{3}1.9477340411.94773404122
5556\frac{5}{6}2.3009758912.30097589144
66112.7182818282.71828182811

Using eqref(52),

S6=h3[f(x0)+4(f(x1)+f(x3)+f(x5))+2(f(x2)+f(x4))+f(x6)]=118[1+4∑i=13e(2i−1)/6+2∑i=12e2i/6+e]=118[1+4(1.181360413+1.648721271+2.300975891)+2(1.395612425+1.947734041)+2.718281828]≈1.718289170.\begin{equation} \begin{aligned} S_6 &=\frac{h}{3}\left[ f(x_0)+4(f(x_1)+f(x_3)+f(x_5)) +2(f(x_2)+f(x_4))+f(x_6) \right] \\ &=\frac{1}{18}\left[ 1+4\sum_{i=1}^{3}e^{(2i-1)/6} +2\sum_{i=1}^{2}e^{2i/6} +e \right] \\ &=\frac{1}{18}\left[ 1+4(1.181360413+1.648721271+2.300975891) +2(1.395612425+1.947734041)+2.718281828 \right] \\ &\approx 1.718289170. \end{aligned} \end{equation}

For f(x)=exf(x)=e^x, we have f(4)(x)=exf^{(4)}(x)=e^x. On [0,1][0,1],

M=max⁡x∈[0,1]∣f(4)(x)∣=e.\begin{equation} M=\max_{x\in[0,1]}\left|f^{(4)}(x)\right|=e. \end{equation}

Using eqref(53),

∣ES∣≤Mnh590=e⋅3⋅(16)590=e233280≈1.16524×10−5.\begin{equation} |E_S|\leq \frac{Mnh^5}{90} =\frac{e\cdot 3\cdot \left(\frac{1}{6}\right)^5}{90} =\frac{e}{233280} \approx 1.16524\times 10^{-5}. \end{equation}

The exact value is e−1e-1, so the true error is

∣(e−1)−S6∣≈∣1.718281828−1.718289170∣=7.34146×10−6.\begin{equation} \left|(e-1)-S_6\right| \approx \left|1.718281828-1.718289170\right| =7.34146\times 10^{-6}. \end{equation}

Splitting [a,b][a, b] into smaller subintervals (making hh smaller) reduces the error for composite rules. For the composite trapezoidal rule,

∣ET∣≤M(b−a)12h2,ET=O(h2)|E_T| \leq \frac{M(b-a)}{12} h^2, \qquad E_T = O(h^2)

Halving hh reduces the error by about a factor of 44. For the composite Simpson’s rule,

∣ES∣≤M(b−a)180h4,ES=O(h4)|E_S| \leq \frac{M(b-a)}{180} h^4, \qquad E_S = O(h^4)

Halving hh reduces the error by about a factor of 1616. Simpson’s rule converges much faster than the trapezoidal rule for smooth functions.

The exponent of hh in the error term is what matters here. Higher order requires (often, but not necessarily) fewer points for the same accuracy (i.e., the function has to be nice relative to the approximation method being used). That’s the whole point, not just getting the answer but getting the right answer while using as little resources as possible (computationally speaking). In practice, we don’t really care (much) about the error term, rather, we care about the order of convergence. The two are pretty much the same thing in a sense (not really, but also, not really still), but the idea is what matters. In practice, there is a lot of nuance and cases which you should check and consider to get a feel for how these work among different functions, and I won’t be going over those cause you already have all the information you need to test it out yourself. But you should have a look and see how these approximation methods are used and in which contexts they are used as well as why some methods aren’t used anymore (generally speaking). \newline \newline Anyway that’s it for now, there are a lot of details that I’ve left out, but that’s what books are for.