Calculus with Code
There is a lot of overlap between the fields of mathematics and programming, and there is computing software that implements mathematics with code, including MATLAB, Mathematica and SageMath.
While such programs are already capable of handling most mathematical tasks, including calculus of course, I'd still like to post some primitive code implementations (i.e. without using thirdparty modules like NumPy and considering floatingpoint errors) of the two major components of calculus, i.e. differentiation and integration, using Python, which I've learnt from my university programming methodology course, here. They also serve as demonstrations of the concept of higherorder functions, which they mathematically are as well.
Here, instead of deriving the symbolic expressions of the derivatives or integrals, the numerical methods of evaluating them are presented.
Differentiation
First Derivative
The following code snippet takes in a (realvalued) function f
and returns a function that takes in a value a
(on which f
is differentiable) and returns the approximated first derivative of f
, applying the first principle of differentiation.
def derivative(f):
h = 10 ** (10) # hardcoded small value for accuracy
return lambda a: (f(a+h)  f(a)) / h
Higherorder Derivatives
To compute the higherorder derivatives, one can use the following code snippet that involves selfcomposition of functions, or more formally: (positive) functional powers.
# Iterative functional composition implementation
def compose(f, pow):
composed = lambda x: x # identity function
for i in range(pow):
composed = (lambda g: lambda x: f(g(x)))(composed)
return composed
def higher_derivative(f, ord):
return lambda x: compose(derivative(f), ord)(x)
degree_five_polynomial = lambda x: x ** 5  x + 1
df, d5f = derivative(degree_five_polynomial), higher_derivative(degree_five_polynomial, 5)
print(df(0)) # 1.000000082740371, actual: 1
print(d5f(0)) # 0.0, actual: 0.0
NewtonRaphson Method
The NewtonRaphson method is an iterative algorithm that makes use of derivatives to approximate a root of a (realvalued) function $f$ (which is differentiable near the root), i.e. the value $a$ for which $f(a)=0$.
$\text{Start from an initial guess } x_0,\; \text{ then for } n \geqslant 0:\; x_{n+1}=x_{n}{\frac {f(x_{n})}{f'(x_{n})}}.$The following code implementation takes in a function f
, an inital guess of the root value init
, and optionally the number of decimal places d
that the estimate should be accurate up to, then outputs the resulting approximated value.
def newtons_method(f, init, *d):
if len(d) > 1:
raise TypeError("Too many arguments")
df = derivative(f) # reusing the derivative function above
def improve(x):
return x(f(x)/df(x))
tolerance = 10 ** (15) # hardcoded small value for accuracy
guess = init
while abs(f(guess)) > tolerance:
guess = improve(guess)
return round(guess, d[0]) if d else guess
square = lambda x: x*x
sqrt_newton = lambda a: newtons_method(lambda t: square(t)a, a/2)
print(sqrt_newton(2)) # 1.4142135623730951
Integration  Definite Integrals: Area under a Curve
Unless otherwise specified, the Riemann version of integrals is assumed to be default here.
The following code snippets approximate the area under the graph of a real function f
between two given values a
and b
, applying the concept of Riemann sum) (rectangle rule).
The midpoint rule, trapezium rule, Simpson's rule and Boole's rule^{1} are covered here.
Midpoint Rule
The midpoint rule approximates the definite integral by evauluating and summing over the values at the midpoints of a number of subintervals that partition the two given endpoints of the integral. It is 3rdorder accurate, i.e. it has an error of order $O(h^3)$.
$\int_{a}^{b}f(x)dx \approx \Delta x\left[f\left(a+{\tfrac {\Delta x}{2}}\right)+f\left(a+{\tfrac {3\Delta x}{2}}\right)+\dots +f\left(b{\tfrac {\Delta x}{2}}\right)\right]$The following code snippet makes use of the sum
function that simulates the behaviour of summation.
def sum(term, a, next, b):
result, curr_idx = 0, a
while curr_idx <= b:
result += term(curr_idx)
curr_idx = next(curr_idx)
return result
def integral_midpoint(f, a, b, dx):
'''
This function takes in the following parameters:
f: function to be integrated,
a: lower limit of integration,
b: upper limit of integration,
dx: the length of each subinterval in the equal partition of the integration interval,
and outputs the value of the integral approximated using the midpoint rule.
'''
add_dx = lambda x: x + dx
return dx * sum(f, a+(dx/2), add_dx, b)
cube = lambda x: x*x*x
print(integral_midpoint(cube, 0, 1, 0.01)) # 0.24998750000000042, exact value: 0.25
Trapezium Rule
The trapezium rule is a technique that approximates the definite integral via linear functions passing through the endpoints of the subintervals, thereby resulting in trapeziums whose areas we can then sum over. It is 3rdorder accurate, i.e. it has an error of order $O(h^3)$.
$\int _{a}^{b}f(x)\,dx \approx {\frac {h}{2}}{\Biggl (}f(x_{0})+2f(x_{1})+2f(x_{2})+2f(x_{3})+\dotsb +2f(x_{n1})+f(x_{n}){\Biggr )} \text{ where } h=\frac{ba}{n}$The sum
function defined as above is reused in the following implementation.
def integral_trapezium(f, a, b, n):
'''
This function takes in the following parameters:
f: function to be integrated,
a: lower limit of integration,
b: upper limit of integration,
n: number of subintervals in the partition of the integration interval,
and outputs the value of the integral approximated using the trapezium rule.
'''
h = (ba)/n
def term(i):
t, c = a + i * h, None
if i == 0 or i == n:
c = 1
else:
c = 2
return c * f(t)
return (h/2) * sum(term, 0, lambda i: i + 1, n)
cube = lambda x: x*x*x
print(integral_trapezium(cube, 0, 1, 50)) # 0.25010000000000004, exact value: 0.25 (set n = 50)
Simpson's Rule
The Simpson's rule is a refinement of the midpoint and trapezium rules in the sense that it has a better order of accuracy; it is 5thorder accurate, i.e. an error of order $O(h^5)$. This rule approximates the integral via a bunch of quadratic functions (parabolas).
The sum
function defined as above is reused (again) in the following implementation.^{2}
def integral_simpson(f, a, b):
'''
This function takes in the following parameters:
f: function to be integrated,
a: lower limit of integration,
b: upper limit of integration,
and outputs the value of the integral approximated using the Simpson rule.
'''
return integral_comp_simpson(f, a, b, 2) # check the "Composite Simpson's Rule" section for the definition of integral_comp_simpson
cube = lambda x: x*x*x
print(integral_simpson(cube, 0, 1)) # 0.25, exact value: 0.25
Composite Simpson's Rule
The composite Simpson's rule is a generalisation of the Simpson's rule that is useful in approximating the definite integral of a function over an interval on which it is not smooth.
One obtains the basic Simpson's rule by setting $n = 2$.
$\int_{a}^{b}f(x)dx\approx{\frac {h}{3}}{\big [}f(x_{0})+4f(x_{1})+2f(x_{2})+4f(x_{3})+2f(x_{4})+\dots +2f(x_{n2})+4f(x_{n1})+f(x_{n}){\big ]} \text{ where } h=\frac{ba}{n}$This was featured as a tutorial question of the aforementioned programming methodology course.
The following code implementation is my solution to the question, which also makes use of the sum
function.
def integral_comp_simpson(f, a, b, n):
'''
This function takes in the following parameters:
f: function to be integrated,
a: lower limit of integration,
b: upper limit of integration,
n: number of subintervals in the partition of the integration interval,
and outputs the value of the integral approximated using the composite Simpson rule.
'''
h = (b  a)/n
def term(i):
t, c = a + i * h, None
if i == 0 or i == n:
c = 1
elif i % 2 == 1:
c = 4
elif i % 2 == 0:
c = 2
return c * f(t)
return (h/3) * sum(term, 0, lambda i: i + 1, n) # sum function defined as above
cube = lambda x: x*x*x
print(integral_comp_simpson(cube, 0, 1, 10)) # 0.25, exact value: 0.25 (set n = 10)
Boole's Rule
The Boole's rule is an improvement of the Simpson's rule with an even higher order of accuracy. It is 7thorder accurate, i.e. an error of order $O(h^7)$. Here, the composite Boole's rule will be used.
This rule approximates the integral via quartic interpolation, using a bunch of quartic functions that passes through five evenly spaced points in each subinterval. Such application of a higherdegree function in turn contributes to the enhanced accuracy of this method.
$\int_{a}^{b}f(x)dx \approx {\frac {2h}{45}}\left(7(f(a)+f(b))+32\left(\sum _{i\in \{1,3,5,\ldots ,N1\}}f(x_{i})\right)+12\left(\sum _{i\in \{2,6,10,\ldots ,N2\}}f(x_{i})\right)+14\left(\sum _{i\in \{4,8,12,\ldots ,N4\}}f(x_{i})\right)\right) \text{ where } h=\frac{ba}{n}$def integral_comp_boole(f, a, b, n):
'''
This function takes in the following parameters:
f: function to be integrated,
a: lower limit of integration,
b: upper limit of integration,
n: number of subintervals in the partition of the integration interval, where n is divisible by 4,
and outputs the value of the integral approximated using the composite Boole's rule.
'''
if n % 4 != 0:
raise ValueError("n must be divisible by 4")
h = (b  a)/n
def term(i):
t, c = a + i * h, None
if i == 0 or i == n:
c = 7
elif i % 4 == 1 or i % 4 == 3:
c = 32
elif i % 4 == 2:
c = 12
elif i % 4 == 0:
c = 14
return c * f(t)
return (2*h/45) * sum(term, 0, lambda i: i + 1, n) # sum function defined as above
cube = lambda x: x*x*x
print(integral_comp_boole(cube, 0, 1, 4)) # 0.25, exact value: 0.25 (set n = 4)
References, further reading
 Chang, X.W. (2023). Numerical Integration. (McGill University COMP 350 Numerical Computing, Fall 2023 Teaching Slide)
 Pflueger, N. (2011). Math 1B, lecture 4: Error bounds for numerical methods.
 Composite Rules. Harvey Mudd College.
Footnotes

The trapezium rule, Simpson's rule and Boole's rule are part of a bigger family of numerical integration methods called NewtonColes formulae. ↩

The fact that the
sum
function can be reused multiple times like this demonstrates the benefit of functional abstraction and higherorder functions in terms of code clarity. ↩