Riemann Integrability, the Fundamental Theorem of Calculus & Adaptive Simpson's Quadrature

While differential calculus analyzes instantaneous rates of change, integral calculus investigates continuous accumulation, total work, mass distribution, and geometric areas. The transition from discrete sums to continuous integrals forms the bedrock of probability theory, physical mechanics, and numerical quadrature algorithms.

1. Rigorous Construction of the Riemann Integral#

The intuitive notion of the area under a continuous curve requires rigorous mathematical formalization. For non-trivial or rapidly oscillating functions, naive definitions of area collapse into paradoxes. The modern formulation of integration rests upon the work of Bernhard Riemann and Gaston Darboux.

Partitions and Subinterval Meshes#

Let f: [a, b] → R be a bounded real-valued function defined on a closed, bounded interval [a, b].

A partition P of [a, b] is a finite, ordered sequence of points:

P = { x0, x1, x2, ..., xn }, where a = x0 < x1 < x2 < ... < xn = b.

The partition divides the interval [a, b] into n contiguous subintervals:

[x0, x1], [x1, x2], ..., [x_{i-1}, x_i], ..., [x_{n-1}, x_n].

The width of the i-th subinterval is:

Δx_i = x_i - x_{i-1}.

The mesh (or norm) of the partition, denoted ||P||, is the width of the largest subinterval:

||P|| = max { Δx_1, Δx_2, ..., Δx_n }.

Lower and Upper Darboux Sums#

Because function f is bounded on [a, b], it attains a well-defined infimum (greatest lower bound) and supremum (least upper bound) on each subinterval [x_{i-1}, x_i]:

m_i = inf { f(x) : x in [x_{i-1}, x_i] }

M_i = sup { f(x) : x in [x_{i-1}, x_i] }

We construct two step-function approximations:

1. The Lower Darboux Sum L(f, P):

L(f, P) = sum(i=1 to n, m_i · Δx_i).

This represents the area of inscribed rectangles lying entirely beneath the curve.

1. The Upper Darboux Sum U(f, P):

U(f, P) = sum(i=1 to n, M_i · Δx_i).

This represents the area of circumscribed rectangles that completely envelop the curve.

Because m_i <= M_i for every subinterval, it is guaranteed that for any partition P:

L(f, P) <= U(f, P).

+-----------------------------------------------------------------------------+
|                  DARBOUX STEP-FUNCTION SUM APPROXIMATION                    |
+-----------------------------------------------------------------------------+
|    y                                                                        |
|    ^                      f(x)                                              |
|    |                     .------.                                           |
|    |             .------'        `---.                                      |
|    |     M_i +---|                   |---+ <--- Upper Bound M_i             |
|    |         |   |                   |   |      (Circumscribed Rectangle)   |
|    |     m_i +...|                   |...+ <--- Lower Bound m_i             |
|    |         |   |                   |   |      (Inscribed Rectangle)       |
|    |         |   |                   |   |                                  |
|    +---------+---+-------------------+---+---------------------------> x     |
|              a  x_{i-1}             x_i  b                                  |
|                 <------ dx_i ------->                                       |
+-----------------------------------------------------------------------------+

The Refinement Lemma#

A partition P is called a refinement of partition P if P ⊆ P (meaning P* contains all the partition points of P plus additional subdivision points).

Theorem (Monotonic Refinement):

If P* is a refinement of P, then:

L(f, P) <= L(f, P) <= U(f, P) <= U(f, P).

Proof Intuition:

Subdividing an existing interval [x_{i-1}, x_i] by introducing a new point x* allows the infimum of each half to rise (or stay equal), increasing the lower sum. Conversely, the supremum of each half falls (or stays equal), decreasing the upper sum. Refining partitions monotonically squeezes the upper and lower sums together!

Corollary: For ANY two arbitrary partitions P1 and P2 of [a, b]:

L(f, P1) <= U(f, P2).

No lower sum can ever exceed any upper sum.

The Riemann Integrability Criteria#

Because the set of all lower sums is bounded above, it possesses a unique supremum. Because the set of all upper sums is bounded below, it possesses a unique infimum:

• Lower Darboux Integral: L_integral = sup { L(f, P) : P is a partition of [a, b] }.
• Upper Darboux Integral: U_integral = inf { U(f, P) : P is a partition of [a, b] }.

Definition: A bounded function f is Riemann integrable on [a, b] if and only if the upper and lower Darboux integrals coincide:

L_integral = U_integral = integral_a^b f(x) dx.

Cauchy's Integrability Criterion#

Theorem: A bounded function f is Riemann integrable on [a, b] if and only if:

For every ε > 0, there exists a partition P of [a, b] such that:

U(f, P) - L(f, P) < ε.

Lebesgue's Criterion for Riemann Integrability#

Theorem (Henri Lebesgue):

A bounded function f on a compact interval [a, b] is Riemann integrable if and only if its set of discontinuities has Lebesgue measure zero.

Classes of Guaranteed Riemann Integrable Functions:

1. Every continuous function on [a, b] is Riemann integrable.
2. Every monotonic (strictly increasing or decreasing) function on [a, b] is Riemann integrable, even if it possesses infinitely many jump discontinuities!
3. Piecewise continuous functions with a finite number of jump discontinuities.

Pathological Non-Integrable Counterexample (The Dirichlet Function):

Define D(x) on [0, 1] by:

D(x) = 1 if x is rational, and 0 if x is irrational.

On any subinterval [x_{i-1}, x_i], no matter how small, there exist both rational and irrational numbers.

Therefore: m_i = 0 and M_i = 1 everywhere.

For every possible partition P:

L(D, P) = sum(0 · Δx_i) = 0 => L_integral = 0.

U(D, P) = sum(1 · Δx_i) = 1 => U_integral = 1.

Because L_integral ≠ U_integral (0 ≠ 1), the Dirichlet function is NOT Riemann integrable. (It requires Lebesgue measure integration).

2. The Fundamental Theorem of Calculus (FTC)#

The Fundamental Theorem of Calculus is one of the greatest intellectual achievements in science, uniting two completely disparate operations—differential rates of change and integral accumulation—into mutually reciprocal operations.

Part 1: Differentiation of Accumulation Functions#

Let f: [a, b] → R be continuous on [a, b]. Define the accumulation function F: [a, b] → R by:

F(x) = integral_a^x f(t) dt.

Theorem (FTC Part 1):

F is continuous on [a, b], differentiable on the open interval (a, b), and its derivative is precisely the original integrand:

F'(x) = d/dx [ integral_a^x f(t) dt ] = f(x).

Analytical Proof:

Evaluate the derivative using the difference quotient:

F'(x) = lim_{h → 0} [ F(x + h) - F(x) ] / h.

Express the numerator as integrals:

F(x + h) - F(x) = integral_a^{x + h} f(t) dt - integral_a^x f(t) dt = integral_x^{x + h} f(t) dt.

Therefore:

F'(x) = lim_{h → 0} (1 / h) · integral_x^{x + h} f(t) dt.

By the Mean Value Theorem for Integrals, because f is continuous on [x, x + h], there exists a point c_h strictly between x and x + h such that:

integral_x^{x + h} f(t) dt = f(c_h) · h.

Substitute this into the difference quotient:

(1 / h) · [ f(c_h) · h ] = f(c_h).

Now take the limit as h → 0. As h approaches zero, the window [x, x + h] pinches down, forcing c_h → x.

Because f is continuous:

lim_{h → 0} f(c_h) = f(x).

Hence, F'(x) = f(x). Differentiation undoes integration!

Part 2: Evaluation of Definite Integrals#

Theorem (FTC Part 2):

If f: [a, b] → R is continuous and G is ANY antiderivative of f on [a, b] (meaning G'(x) = f(x) for all x in (a, b)), then:

integral_a^b f(x) dx = G(b) - G(a) = [ G(x) ]_a^b.

Proof:

Let F(x) = integral_a^x f(t) dt. By Part 1, F'(x) = f(x).

Because G'(x) = f(x) as well, F(x) and G(x) share the same derivative.

By the Mean Value Theorem, two functions with identical derivatives differ by a constant:

F(x) = G(x) + C.

At x = a:

F(a) = integral_a^a f(t) dt = 0 => G(a) + C = 0 => C = -G(a).

Therefore:

F(x) = G(x) - G(a).

Evaluating at x = b:

F(b) = integral_a^b f(t) dt = G(b) - G(a).

The Leibniz Integral Rule: Differentiation Under the Integral Sign#

When an integral has variable limits and a parameter inside the integrand:

I(x) = integral_{u(x)}^{v(x)} f(x, t) dt.

The Leibniz Rule states:

dI/dx = f(x, v(x)) · v'(x) - f(x, u(x)) · u'(x) + integral_{u(x)}^{v(x)} [ ∂f(x, t) / ∂x ] dt.

This formula combines the chain rule with the Fundamental Theorem of Calculus and is widely used in physics and machine learning loss estimation.

3. Analytical Integration Techniques: Inverting Differential Operators#

Integration techniques are systematic methods for reversing standard differential rules.

1. Integration by Substitution (u-Substitution)#

Substitution inverts the Chain Rule:

Since d/dx [ F(g(x)) ] = f(g(x)) · g'(x), integrating both sides yields:

integral f(g(x)) · g'(x) dx = integral f(u) du, where u = g(x) and du = g'(x) dx.

For Definite Integrals:

integral_a^b f(g(x)) · g'(x) dx = integral_{g(a)}^{g(b)} f(u) du.

The integration limits MUST be transformed from x-coordinates to u-coordinates!

2. Integration by Parts#

Integration by parts inverts the Product Rule:

d/dx [ u(x) · v(x) ] = u'(x) · v(x) + u(x) · v'(x).

Rearranging and integrating:

integral u · dv = u · v - integral v · du.

For Definite Integrals:

integral_a^b u(x) · v'(x) dx = [ u(x) · v(x) ]_a^b - integral_a^b v(x) · u'(x) dx.

The LIATE Priority Rule for Selecting u(x):

To select which factor to set as u(x) (which will be differentiated) versus dv (which will be integrated):

1. L: Logarithmic functions (ln(x))
2. I: Inverse trigonometric functions (arctan(x), arcsin(x))
3. A: Algebraic polynomials (x^2, x^3)
4. T: Trigonometric functions (sin(x), cos(x))
5. E: Exponential functions (e^x, 2^x)

3. Partial Fraction Decomposition#

Used for rational functions P(x) / Q(x) where degree(P) < degree(Q):

• Distinct linear factors: A / (x - r) + B / (x - s).
• Repeated linear factors: A / (x - r) + B / (x - r)^2 + ...
• Irreducible quadratic factors: (A x + B) / (x^2 + p x + q).

4. Numerical Quadrature: From Riemann Sums to Adaptive Simpson's Rule#

Many functions possess antiderivatives that cannot be expressed in terms of elementary functions (e.g., e^(-x^2), sin(x)/x, sqrt(1 - k^2 sin^2(x))). Evaluating definite integrals of such functions requires numerical quadrature.

The Newton-Cotes Quadrature Family#

Newton-Cotes methods approximate f(x) by dividing [a, b] into n subintervals of uniform step size h = (b - a) / n and fitting interpolating polynomials:

1. Left and Right Riemann Sums:

Fit piecewise constants. Error = O(h).

1. Trapezoidal Rule:

Fits piecewise straight lines across each subinterval:

integral_a^b f(x) dx ≈ (h / 2) · [ f(x0) + 2 f(x1) + 2 f(x2) + ... + 2 f(x_{n-1}) + f(xn) ].

Global Truncation Error:

E_T = -[ (b - a)^3 / (12 · n^2) ] · f''(ξ) = O(h^2).

1. Composite Simpson's 1/3 Rule:

Fits piecewise quadratic parabolas across adjacent pairs of subintervals (requiring an even number of subdivisions n):

integral_a^b f(x) dx ≈ (h / 3) · [ f(x0) + 4 sum_{i odd} f(xi) + 2 sum_{i even} f(xi) + f(xn) ].

Global Truncation Error of Simpson's Rule:

E_S = -[ (b - a)^5 / (180 · n^4) ] · f^(4)(ξ) = O(h^4)

for some point ξ in [a, b].

Why Simpson's Rule is Unusually Powerful:

Even though Simpson's Rule fits second-degree parabolas, its error depends on the FOURTH derivative f^(4)(ξ)!

Because the fourth derivative of any cubic polynomial ax^3 + bx^2 + cx + d is identically zero:

Simpson's Rule integrates all cubic polynomials exactly with zero truncation error!

Adaptive Simpson's Quadrature#

In scientific computing, functions often have regions with steep spikes or rapid oscillations alongside smooth, flat plains. Using a uniform small step size h everywhere wastes computational cycles.

Adaptive Quadrature dynamically refines the step size:

1. Compute Simpson's approximation over the whole interval [a, b]: S_whole.
2. Bisect at midpoint c = (a + b) / 2. Compute Simpson's approximations on left [a, c] (S_left) and right [c, b] (S_right).
3. By Richardson extrapolation, the local error estimate is:

Error_est = |S_left + S_right - S_whole| / 15.

1. If Error_est < Tolerance: Return (S_left + S_right + Error_est).
2. If Error_est >= Tolerance: Recursively call Adaptive Simpson on [a, c] and [c, b] with halved tolerance.

+-----------------------------------------------------------------------------+
|                     ADAPTIVE QUADRATURE TREE REFINEMENT                     |
+-----------------------------------------------------------------------------+
|                                  [ a, b ]                                   |
|                        (Check Error_est < Tolerance)                        |
|                                 /            \                              |
|                         Fails: Subdivide      Passes: Accept                |
|                               /                \                            |
|                       [ a, c ]                  [ c, b ]                    |
|                (High Curvature Region)    (Flat Smooth Region)              |
|                     /          \           Accepted Directly!               |
|                 Subdivide    Subdivide                                      |
+-----------------------------------------------------------------------------+

5. Scientific Computing & Physics Engine Applications#

Integration is the fundamental computational engine for physical simulation, probability, and signal processing.

1. Physics Engines: Symplectic Verlet Integration#

In video game physics engines (Unreal Engine Chaos, PhysX) and molecular dynamics (GROMACS, LAMMPS), simulating Newton's second law F = m · (d^2 x / dt^2) requires integrating accelerations into velocities and positions.

• Naive Explicit Euler: x_{k+1} = x_k + v_k · Δt, v_{k+1} = v_k + a_k · Δt.

Explicit Euler fails to conserve mechanical energy. In an orbiting planet or spring-mass system, numerical error adds artificial energy at every timestep, causing orbits to spiral outward and systems to explode.

• Velocity Verlet (Symplectic Integrator):

x_{k+1} = x_k + v_k · Δt + (1/2) · a_k · Δt^2

v_{k+1} = v_k + (1/2) · (a_k + a_{k+1}) · Δt

Verlet integration is symplectic: it preserves phase space volume and conserves total energy over millions of simulation timesteps without drift.

2. Probability & Continuous Random Variables: CDF Accumulation#

In statistical machine learning, continuous random variables are defined by their Probability Density Function (PDF) p(x). The probability that a random variable X falls in range [a, b] is given by the definite integral:

P(a <= X <= b) = integral_a^b p(x) dx.

The Cumulative Distribution Function (CDF) is the accumulation function:

F(x) = integral_{-infinity}^x p(t) dt.

By the Fundamental Theorem of Calculus, the PDF is the exact derivative of the CDF:

p(x) = d/dx [ F(x) ].

For the Gaussian Normal Distribution, the integral cannot be evaluated in elementary terms and defines the standard Error Function:

erf(x) = (2 / sqrt(π)) · integral_0^x e^(-t^2) dt.

3. Digital Audio Processing: Continuous Convolution Integrals#

In audio DSP and telecommunications, linear time-invariant (LTI) systems filter input signals x(t) through an impulse response h(t) via the continuous convolution integral:

y(t) = (x * h)(t) = integral_{-infinity}^infinity x(τ) · h(t - τ) dτ.

This integral computes the running overlap of the incoming signal with the system's acoustic response, modeling realistic reverb, equalization, and echo cancellation in real time.

6. Common Architectural Pitfalls & Misconceptions#

• Pitfall 1: Blindly Integrating Across Vertical Asymptotes (Improper Integrals).

Evaluating integral_{-1}^1 (1 / x^2) dx by naive antiderivative:

G(x) = -1/x. Then G(1) - G(-1) = (-1/1) - (-1/(-1)) = -1 - 1 = -2.

This answer is completely wrong and physically absurd: 1/x^2 is positive everywhere, so its integral cannot be negative! The integrand diverges to infinity at x = 0; the integral is improper and diverges to +∞.

• Pitfall 2: Forgetting to Transform Integration Limits in u-Substitution.

When evaluating integral_0^1 x · sqrt(1 - x^2) dx with u = 1 - x^2:

If you forget to change the limits from x in [0, 1] to u in [1, 0], and evaluate [-(1/3) u^(3/2)] from 0 to 1, you will get the opposite sign. Always transform the limits immediately upon substitution.

• Pitfall 3: High-Order Newton-Cotes Quadrature Instability (Runge's Phenomenon).

Novice engineers often assume that fitting an 8th-degree or 10th-degree polynomial across equally spaced points will produce extraordinary accuracy. In reality, high-order Newton-Cotes formulas have large negative weights that cause catastrophic roundoff oscillations. Piecewise low-order methods (Composite Simpson or Gauss-Legendre quadrature) are vastly superior.

• Pitfall 4: Neglecting the Trailing Minus Sign in Integration by Parts.

The integration by parts formula is integral u dv = u · v - integral v du. When v du contains a negative sign (e.g., differentiating cos(x) to -sin(x)), the double negative becomes positive. Dropping this sign is the most frequent algebraic source of error.

7. Summary & Checkpoint Exercises#

Solidify your comprehension of Riemann sums, the Fundamental Theorem of Calculus, integration techniques, and numerical quadrature with these four checkpoint problems. Full step-by-step mathematical solutions are provided.

Problem 1 (Calculating Upper and Lower Darboux Sums):

Consider the function f(x) = x^2 on the interval [0, 2].

Divide [0, 2] into n subintervals of uniform width Δx = 2 / n.

1. Formulate the Lower Darboux Sum L(f, P_n) and Upper Darboux Sum U(f, P_n).
2. Evaluate their limits as n → infinity using the summation formula sum(i=1 to n, i^2) = [ n(n + 1)(2n + 1) ] / 6.
3. Prove that f(x) is Riemann integrable on [0, 2] and find the exact value of the integral.

Solution:

Part 1: Formulate Darboux Sums.

Interval [0, 2] partitioned into n subintervals of width Δx = 2/n.

Partition points: x_i = i · (2/n) for i = 0, 1, 2, ..., n.

Subintervals: [ x_{i-1}, x_i ] = [ (i - 1)(2/n), i(2/n) ].

Because f(x) = x^2 is strictly increasing on [0, 2]:

• Infimum on i-th subinterval occurs at the left endpoint:

m_i = f(x_{i-1}) = [ (i - 1)(2/n) ]^2 = (4 / n^2) · (i - 1)^2.

• Supremum on i-th subinterval occurs at the right endpoint:

M_i = f(x_i) = [ i(2/n) ]^2 = (4 / n^2) · i^2.

Assemble the Sums:

L(f, P_n) = sum(i=1 to n, m_i · Δx) = sum(i=1 to n, [ (4 / n^2) · (i - 1)^2 ] · (2 / n))

= (8 / n^3) · sum(i=1 to n, (i - 1)^2) = (8 / n^3) · sum(k=0 to n-1, k^2).

U(f, P_n) = sum(i=1 to n, M_i · Δx) = sum(i=1 to n, [ (4 / n^2) · i^2 ] · (2 / n))

= (8 / n^3) · sum(i=1 to n, i^2).

Part 2: Evaluate Limits as n → infinity.

Using the sum of squares formula:

• Upper Sum:

sum(i=1 to n, i^2) = [ n(n + 1)(2n + 1) ] / 6 = [ 2n^3 + 3n^2 + n ] / 6.

U(f, P_n) = (8 / n^3) · [ (2n^3 + 3n^2 + n) / 6 ] = (8/6) · [ 2 + 3/n + 1/n^2 ] = (4/3) · [ 2 + 3/n + 1/n^2 ].

Taking the limit as n → infinity:

lim_{n → infinity} U(f, P_n) = (4/3) · (2 + 0 + 0) = 8/3.

• Lower Sum:

sum(k=0 to n-1, k^2) = [ (n - 1)(n)(2(n - 1) + 1) ] / 6 = [ (n - 1)(n)(2n - 1) ] / 6 = [ 2n^3 - 3n^2 + n ] / 6.

L(f, P_n) = (8 / n^3) · [ (2n^3 - 3n^2 + n) / 6 ] = (4/3) · [ 2 - 3/n + 1/n^2 ].

Taking the limit as n → infinity:

lim_{n → infinity} L(f, P_n) = (4/3) · (2 - 0 + 0) = 8/3.

Part 3: Riemann Integrability and Integral Value.

Because lim_{n → infinity} L(f, P_n) = lim_{n → infinity} U(f, P_n) = 8/3:

The upper and lower Darboux integrals are identical:

sup_P L(f, P) = inf_P U(f, P) = 8/3.

Therefore, f(x) = x^2 is strictly Riemann integrable on [0, 2], and:

integral_0^2 x^2 dx = 8/3.

(Verification via FTC: [ x^3 / 3 ]_0^2 = 2^3 / 3 - 0 = 8/3).

Problem 2 (Leibniz Rule for Differentiating an Integral with Variable Limits):

Find the derivative dF/dx of the function:

F(x) = integral_{x^2}^{x^3} e^(-t^2) dt.

Solution:

Step 1: Identify Functions for Leibniz Integral Rule.

Recall the general Leibniz rule:

d/dx [ integral_{u(x)}^{v(x)} f(x, t) dt ] = f(x, v(x)) · v'(x) - f(x, u(x)) · u'(x) + integral_{u(x)}^{v(x)} [ ∂f/∂x ] dt.

Here:

• Integrand: f(x, t) = e^(-t^2). Notice f does not depend explicitly on x, so ∂f/∂x = 0.
• Upper limit: v(x) = x^3 => v'(x) = 3 x^2.
• Lower limit: u(x) = x^2 => u'(x) = 2 x.

Step 2: Apply the Formula.

1. Evaluate integrand at upper limit:

f(v(x)) = e^(-(x^3)^2) = e^(-x^6).

Multiply by v'(x):

e^(-x^6) · (3 x^2) = 3 x^2 · e^(-x^6).

1. Evaluate integrand at lower limit:

f(u(x)) = e^(-(x^2)^2) = e^(-x^4).

Multiply by u'(x):

e^(-x^4) · (2 x) = 2 x · e^(-x^4).

1. Interior partial derivative integral:

integral_{x^2}^{x^3} 0 dt = 0.

Step 3: Combine Terms.

dF/dx = 3 x^2 · e^(-x^6) - 2 x · e^(-x^4).

This is the exact analytical derivative.

Problem 3 (Repeated Integration by Parts for a Polynomial-Exponential Function):

Evaluate the definite integral:

I = integral_0^1 x^2 · e^(2x) dx.

Solution:

Step 1: First Integration by Parts.

Using LIATE, choose:

u = x^2 => du = 2x dx

dv = e^(2x) dx => v = (1/2) · e^(2x)

Apply formula integral u dv = u · v - integral v du:

integral x^2 · e^(2x) dx = (1/2) x^2 · e^(2x) - integral (1/2) e^(2x) · (2x) dx

= (1/2) x^2 · e^(2x) - integral x · e^(2x) dx.

Step 2: Second Integration by Parts on integral x · e^(2x) dx.

Choose:

u = x => du = dx

dv = e^(2x) dx => v = (1/2) · e^(2x)

Apply formula:

integral x · e^(2x) dx = (1/2) x · e^(2x) - integral (1/2) e^(2x) dx

= (1/2) x · e^(2x) - (1/4) e^(2x).

Step 3: Combine the Antiderivative.

Substitute back:

G(x) = (1/2) x^2 · e^(2x) - [ (1/2) x · e^(2x) - (1/4) e^(2x) ]

= [ (1/2) x^2 - (1/2) x + 1/4 ] · e^(2x).

Step 4: Evaluate Definite Limits from 0 to 1.

• At upper limit x = 1:

G(1) = [ (1/2)(1)^2 - (1/2)(1) + 1/4 ] · e^(2(1))

= [ 1/2 - 1/2 + 1/4 ] · e^2 = (1/4) · e^2.

• At lower limit x = 0:

G(0) = [ 0 - 0 + 1/4 ] · e^(0) = (1/4) · 1 = 1/4.

Definite Integral Value:

I = G(1) - G(0) = (1/4) · e^2 - 1/4 = (e^2 - 1) / 4.

Evaluating numerically: (7.389056 - 1) / 4 = 6.389056 / 4 ≈ 1.597264.

Problem 4 (Composite Simpson's Rule and Error Bound Calculation):

Consider approximating the integral:

I = integral_0^2 (x^4 + 1) dx

using Composite Simpson's 1/3 Rule with n = 4 subintervals.

1. Calculate the step size h and the partition points x0, x1, x2, x3, x4.
2. Compute the Simpson's Rule approximation S4.
3. Compute the exact integral via FTC and determine the actual error.
4. Verify the theoretical error bound formula E_S = -[ (b - a)^5 / (180 · n^4) ] · f^(4)(ξ).

Solution:

Part 1: Setup Partition.

Interval [a, b] = [0, 2], n = 4.

Step size: h = (b - a) / n = (2 - 0) / 4 = 0.5.

Partition points:

x0 = 0.0, x1 = 0.5, x2 = 1.0, x3 = 1.5, x4 = 2.0.

Evaluate f(x) = x^4 + 1 at partition points:

• f(x0) = f(0.0) = 0^4 + 1 = 1.0
• f(x1) = f(0.5) = (0.5)^4 + 1 = 0.0625 + 1 = 1.0625
• f(x2) = f(1.0) = 1^4 + 1 = 1 + 1 = 2.0
• f(x3) = f(1.5) = (1.5)^4 + 1 = 5.0625 + 1 = 6.0625
• f(x4) = f(2.0) = 2^4 + 1 = 16 + 1 = 17.0

Part 2: Compute Simpson's Approximation S4.

Simpson's 1/3 Formula:

S4 = (h / 3) · [ f(x0) + 4 f(x1) + 2 f(x2) + 4 f(x3) + f(x4) ]

S4 = (0.5 / 3) · [ 1.0 + 4(1.0625) + 2(2.0) + 4(6.0625) + 17.0 ]

= (1 / 6) · [ 1.0 + 4.25 + 4.0 + 24.25 + 17.0 ]

= (1 / 6) · [ 50.5 ]

= 50.5 / 6 ≈ 8.416667.

Part 3: Compute Exact Integral and Actual Error.

Antiderivative:

G(x) = x^5 / 5 + x.

G(2) - G(0) = (2^5 / 5 + 2) - 0 = 32/5 + 2 = 6.4 + 2 = 8.4.

Exact Integral: I_exact = 8.4 = 42/5.

Actual Error:

E_actual = I_exact - S4 = 8.4 - 8.416667 = -0.016667 = -1/60.

Part 4: Verify Theoretical Error Bound.

Compute derivatives of f(x) = x^4 + 1:

• f'(x) = 4 x^3
• f''(x) = 12 x^2
• f'''(x) = 24 x
• f^(4)(x) = 24 (constant for all x!).

Because f^(4)(x) = 24 is constant, f^(4)(ξ) = 24 for any point ξ in [0, 2].

Theoretical error:

E_S = -[ (b - a)^5 / (180 · n^4) ] · f^(4)(ξ)

= -[ (2 - 0)^5 / (180 · 4^4) ] · 24

= -[ 32 / (180 · 256) ] · 24

= -[ 32 · 24 ] / [ 46080 ]

= -768 / 46080.

Divide numerator and denominator by 768:

46080 / 768 = 60.

E_S = -1/60 ≈ -0.016667.

The theoretical error formula matches the actual error to the exact decimal!