SymPy Example: Deriving Simpson’s Rule

from sympy import init_session
init_session(use_latex="mathjax")

IPython console for SymPy 1.14.0 (Python 3.14.3-64-bit) (ground types: python)

These commands were executed:
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at https://docs.sympy.org/1.14.0/

We want to fit a parabola to 3 points: \((x_0, f_0)\), \((x_1, f_1)\), \((x_2, f_2)\) and then integrate under this fit. We’ll assume that these points are equally spaced.

We’ll choose a polynomial of the form:

\[f(x) = a (x - x_0)^2 + b (x - x_0) + c\]

a, b, c = symbols("a b c")
x0 = symbols("x_0")
dx = symbols(r"\Delta{}x")
f0, f1, f2 = symbols("f_0 f_1 f_2")

f = a * (x - x0)**2 + b * (x - x0) + c

Now we can get 3 constraints by evaluating this at our 3 points

constraint_0 = Eq(f.subs(x, x0), f0)
constraint_0

\[\displaystyle c = f_{0}\]

constraint_1 = Eq(f.subs(x, x0 + dx), f1)
constraint_1

\[\displaystyle \Delta{}x^{2} a + \Delta{}x b + c = f_{1}\]

constraint_2 = Eq(f.subs(x, x0 + 2 * dx), f2)
constraint_2

\[\displaystyle 4 \Delta{}x^{2} a + 2 \Delta{}x b + c = f_{2}\]

This is a linear system that we can solve algebraically

coeffs = solve([constraint_0, constraint_1, constraint_2],
               [a, b, c])
coeffs

\[\displaystyle \left\{ a : \frac{f_{0} - 2 f_{1} + f_{2}}{2 \Delta{}x^{2}}, \ b : \frac{- 3 f_{0} + 4 f_{1} - f_{2}}{2 \Delta{}x}, \ c : f_{0}\right\}\]

This shows that our polynomial is

\[f(x) = \frac{f_0 - 2 f_1 + f_2}{2 \Delta x^2} (x - x_0)^2 + \frac{-3 f_0 + 4 f_1 - f_2}{2 \Delta x} (x - x_0) + f_0\]

Now we can integrate under this from \([x_0, x_2]\) (or equivalently \([x_0, x_0 + 2 \Delta x]\)

I = integrate(f.subs({a: coeffs[a], b: coeffs[b], c: coeffs[c]}),
              [x, x0, x0+2*dx])
simplify(I)

\[\displaystyle \frac{\Delta{}x \left(f_{0} + 4 f_{1} + f_{2}\right)}{3}\]

This shows that the integral is:

\[I = \int_a^b f(x) dx \approx \frac{\Delta x}{3} \left [ f(a) + 4 f\left(\frac{a+b}{2}\right ) + f(b) \right ]\]