Numerical differentiation and integration

State the target and the inputs

Consider a supplied position model with values 0.65 m at 0.5 s and 2.75 m at 1.5 s along one fixed axis. Across that 1 s separation, the rate is (2.75 − 0.65)m/(1 s) = 2.1 m/s. Does it also give the rate at exactly 1 s? The two values alone do not settle that question.

A derivative describes a local rate; a numerical difference uses values at finitely separated inputs. An integral describes signed accumulation over an interval; numerical quadrature replaces it with a weighted sum. The finite result and the calculus target are distinct, even when they happen to agree.

Let t be time in seconds and q(t) a scalar position in metres on one fixed axis. The target q′(t) has units m/s. Separately, a supplied scalar velocity v(t), in m/s, has signed accumulation I = ∫ₐᵇ v(t) dt, in metres, for a < b. This is displacement along the stated axis; distance would use ∫ₐᵇ |v(t)| dt instead.

Values may come from an exactly specified function or from records with errors. A formula does not supply missing information about the function between records. Before calculating, state the permitted domain, nodes, units, target and any bounds on derivatives or input errors.

For preparation, use Derivatives, Taylor series, and local approximation for derivatives and Taylor notation, and Integrals and multivariable calculus for signed integrals. On a first reading, follow the finite calculations and return to the smoothness-based error bounds after those lessons.

Construct a one-sided difference

For a positive time step h, the forward difference is D_F(t,h) = [q(t + h) − q(t)]/h. The backward difference is D_B(t,h) = [q(t) − q(t − h)]/h. Each uses a separation h and has rate units. All required nodes must lie in the supplied domain, and h must be nonzero.

If q has two continuous derivatives on the interval used, Taylor's theorem bounds the forward formula error by |D_F − q′(t)| ≤ M₂h/2, where M₂ bounds |q″| on [t, t + h]. For the backward formula use a bound on [t − h, t]. M₂ has units m/s², so the error bound has units m/s.

These are exact-input, exact-arithmetic formula bounds. They do not include errors in the values or operations. At the beginning of a restricted record interval, a forward construction may be available while a centered construction would need a missing earlier value.

Center the stencil

The centered difference is D_C(t,h) = [q(t + h) − q(t − h)]/(2h). Its endpoints are separated by 2h, not h. With three continuous derivatives on [t − h, t + h], the bound |D_C − q′(t)| ≤ M₃h²/6 holds when M₃ bounds |q‴| there. M₃ has units m/s³.

Taylor expansions on the two sides cancel the equal quadratic terms in the numerator. Dividing by 2h leaves the derivative and a remainder bounded by M₃h²/6. With the same uniform derivative bound, halving h quarters this upper bound; it need not quarter the actual error of an arbitrary function.

For an original exact model on 0 s ≤ t ≤ 2 s, supply q(t) = 0.2 m + (0.8 m/s)t + (0.4 m/s³)t³. At t = 1 s, q′ = 2 m/s and q‴ = 2.4 m/s³. With h = 0.5 s, the forward, backward and centered differences are 2.7, 1.5 and 2.1 m/s. The centered error is +0.1 m/s; at h = 0.25 s it is +0.025 m/s. These values follow from the stipulated polynomial.

The figure compares the joining secant with the tangent at 1 s. The secant uses the two outer positions and has slope 2.1 m/s; the tangent touches the curve at the middle position and has slope 2 m/s. The middle position is 1.40 m, while the secant at that time is 1.70 m: a centered secant need not pass through the center value.

Include errors in supplied values

Suppose the two centered endpoint values supplied to the calculation are q̃(t ± h) = q(t ± h) + ε±, with |ε±| ≤ δ_q. Times are exact, δ_q is a length, and the subsequent arithmetic is exact in this model. No probability distribution, independence or confidence interpretation is supplied.

The additional difference error is (ε+ − ε−)/(2h), whose magnitude is at most δ_q/h. Combining it with the smooth-function bound gives |D̃_C − q′(t)| ≤ M₃h²/6 + δ_q/h. The first term decreases as h shrinks; the second increases when a fixed positive δ_q is retained.

For the polynomial above with δ_q = 0.010 m, steps 0.5, 0.25 and 0.125 s give total bounds 0.120, 0.065 and 0.08625 m/s. The middle step has the smallest bound among these three choices. This is a comparison under the declared bounds, not a universal optimal step or a measurement-uncertainty estimate.

Build a signed trapezoidal sum

For integration, divide [a, b] into n equal panels, with positive integer n, h = (b − a)/n and tⱼ = a + jh for j = 0,…,n. A line segment joining two adjacent velocity values has signed integral h[v(tⱼ) + v(tⱼ₊₁)]/2. Adding panels gives Tₙ = h[v(a)/2 + Σⱼ₌₁ⁿ⁻¹ v(tⱼ) + v(b)/2]. For n = 1 the interior sum is empty.

The interior values receive a full h weight because each belongs to two adjacent panels; each endpoint receives h/2. The weights have time units and the weighted velocity values give lengths. Negative values contribute negatively. Use the actual panel widths for unequal nodes rather than inserting them into this equal-grid formula.

For a separate original velocity model on 0 s ≤ t ≤ 2 s, supply v(t) = 1 m/s − (0.5 m/s³)t². It is not q′ from the earlier position model. Its exact signed integral is I = 2/3 m. With n = 2 and h = 1 s, the node values are 1, 0.5 and −1 m/s, giving T₂ = 0.5 m. A negative final velocity is retained in this displacement calculation.

Bound the quadrature error

If v has two continuous derivatives on [a, b] and K₂ bounds |v″| there, then |Tₙ − I| ≤ (b − a)K₂h²/12. This is an exact-value formula bound. K₂ has units m/s³, so time length times h² times K₂ has length units.

On each panel, the deviation between a function and its joining line can be bounded using the second derivative. Integrating that deviation gives K₂h³/12 per panel; summing n panels gives the stated bound. A function affine in time has zero second derivative and its trapezoidal integral is exact.

For the supplied quadratic velocity, v″ = −1 m/s³ everywhere, so K₂ = 1 m/s³. Every trapezoidal result lies below the exact signed integral. Here the signed error is exactly −(b − a)(1 m/s³)h²/12, so its magnitude equals the bound. The table uses metres throughout its result and error columns.

Check what the values can establish

If supplied node velocities ṽ(tⱼ) = v(tⱼ) + εⱼ obey |εⱼ| ≤ δ_v and arithmetic is exact, the positive trapezoidal weights sum to b − a. Therefore |T̃ₙ − Tₙ| ≤ (b − a)δ_v. A common signed offset can attain this bound; adding more nodes does not automatically remove it.

For δ_v = 0.020 m/s over 2 s, the value-error bound is 0.040 m. At n = 4, combining it with the formula bound gives |T̃₄ − I| ≤ 1/24 m + 0.040 m = 49/600 m, about 0.08167 m. This bound cannot certify a requested 0.050 m tolerance, although a particular realized error might be smaller.

Two endpoint values alone also do not determine the integral. On the same interval, w(t) = v(t) + (0.75 m/s³)t(2 s − t) shares v's endpoint values but has integral 5/3 m instead of 2/3 m. Its one-panel trapezoidal result is still 0 m. A bound on the unseen curvature must come from a stated model or other evidence, not from agreement at the endpoints.

Interpret a refinement check

Refinement compares calculations with more closely spaced nodes. For exact values with the same applicable derivative bounds, these centered and trapezoidal bounds scale as h². If input errors, timing errors or finite arithmetic matter, account for them before interpreting the comparison. A small difference between two computed results is not by itself a bound on either one's error.

The smoothness hypotheses matter. A corner can have no derivative even though a centered quotient returns a finite number. The trapezoidal bound used here needs a continuous second derivative; a singularity or discontinuity requires a different analysis. Stay inside the specified domain and check that the required values and arithmetic are defined.

Record the nodes, method, arithmetic assumptions, target units, applicable bound and tolerance decision with the result. Numerical differentiation and quadrature act on supplied function values. The next subject, numerical solutions of ordinary differential equations, addresses an unknown function constrained by a differential equation and initial or boundary data.

Worked examples

Try it yourself

1. For the supplied cubic position at t = 1 s, compute D_F with h = 0.25 s. Is a centered stencil with this positive h available at t = 0 s using only the stated domain?

   Show the answer and reasoning

   D_F(1 s, 0.25 s) = 2.325 m/s. The centered stencil at 0 s is unavailable because it would require q(−0.25 s).

   q(1.25 s) = 1.98125 m and q(1 s) = 1.40 m, so divide their 0.58125 m difference by 0.25 s. A domain restriction is not permission to extrapolate the missing value.

2. For the same cubic at t = 1 s, use h = 0.125 s to find D_C and its signed exact-value error.

   Show the answer and reasoning

   D_C = 2.00625 m/s and D_C − q′ = +0.00625 m/s.

   For this cubic, D_C − q′ = (0.4 m/s³)h² exactly. Substituting h = 0.125 s gives 0.00625 m/s. Both required nodes remain in [0 s, 2 s].

3. Each of the two values in a forward difference with h = 0.25 s has an additive position error of magnitude at most 0.010 m. With exact times and arithmetic, what bound applies to the difference error caused only by those value errors?

   Show the answer and reasoning

   The value-error contribution is bounded by 2δ_q/h = 0.080 m/s.

   The numerator error is ε₁ − ε₀ and has magnitude at most 0.020 m. Dividing by 0.25 s gives 0.080 m/s. This is separate from the forward formula error relative to the derivative.

4. For the supplied quadratic velocity, give T₂, T₄ and their absolute exact-value errors. Does quartering of these errors establish the same behavior for every supplied data set?

   Show the answer and reasoning

   T₂ = 1/2 m and T₄ = 5/8 m; absolute errors are 1/6 m and 1/24 m. The factor of four is established for this exact quadratic case.

   I = 2/3 m. The constant second derivative makes the trapezoidal error exactly proportional to h² here. For another smooth function the general bound scales with h² under a common K₂; its realized error need not follow that exact ratio, and input errors require separate treatment.

5. At n = 4, also suppose each supplied velocity value differs from the model by at most 0.020 m/s. What total absolute bound follows with exact times and arithmetic, and does it certify 0.050 m?

   Show the answer and reasoning

   The total bound is 1/24 m + 0.040 m = 49/600 m, about 0.08167 m. It does not certify 0.050 m.

   Positive trapezoidal weights sum to 2 s, so node-value errors contribute at most (2 s)(0.020 m/s) = 0.040 m. Add this to the formula bound by the triangle inequality. The failure to certify does not specify the unknown realized error.

6. On [0 s, 2 s], compare v(t) with w(t) = v(t) + (0.75 m/s³)t(2 s − t). Check their endpoint values, one-panel trapezoidal results and exact integrals.

   Show the answer and reasoning

   Both have endpoints 1 and −1 m/s and one-panel result 0 m. Their exact integrals are 2/3 m for v and 5/3 m for w.

   The added term vanishes at both endpoints, yet its integral is (0.75 m/s³)[(2 s)t²/2 − t³/3] evaluated from 0 s to 2 s, which equals 1 m. Agreement at the two nodes therefore does not determine the signed integral; the functions have different curvature.

What to carry forward

• A finite difference or quadrature sum must identify the target quantity, units, nodes and valid domain.
• Smoothness and derivative bounds control exact-value formula error; they are assumptions requiring support.
• Shrinking a difference step can amplify fixed errors in its supplied endpoint values.
• Trapezoidal integration preserves signed contributions; its positive weights bound uniform node errors without assuming independence.
• Report the applicable total bound and tolerance decision, and distinguish failure to certify from an observed failure.