Python Calculus Toolkit

A Complete Code Explanation & Calculus Reference

Covers: Definite Integrals | Numerical Methods | Applications

1. Overview
This Python script is a calculus computation toolkit that solves 10 exercises involving
definite integrals, numerical integration methods, and real-world applications. It is built on
three powerful libraries that together cover both symbolic (exact) and numerical
(approximate) mathematics.

Library Type Purpose

sympy Symbolic math Computes exact integrals, simplifications, and
mathematical expressions — like a CAS
(Computer Algebra System)
numpy Numerical math Fast array operations and approximate
calculations using floating-point arithmetic
matplotlib Visualization Plots graphs of mathematical functions for visual
analysis

The script is structured as a sequence of exercises (Exercise 1 through Exercise 10), each
targeting a specific calculus concept. Before the exercises, several helper functions are
defined to avoid code repetition.
2. Helper Functions
These functions are defined at the top of the file and reused throughout all exercises.
Understanding them is essential before reading the exercises.

2.1 show(name, items)

def show(name, items):
print(f"\n{name}:")
for k, v in items:
print(f"({k}) {[Link](v)}")

Purpose: A pretty-printer for displaying exercise results. It accepts a list of (label,

expression) tuples and prints each one after simplifying it with SymPy.

[Link](v): SymPy's simplification function reduces an expression to its simplest form.

For example, sin(x)^2 + cos(x)^2 simplifies to 1. Without this call, results can look more
complex than needed.

CS Concept: The parameter items is a list of tuples — a common Python pattern for pairing
related data. The for k, v in items line uses tuple unpacking, which is cleaner than using
items[i][0] and items[i][1].

2.2 plot_expr(expr, var, start, end, title)

def plot_expr(expr, var, start, end, title):
f = [Link](var, expr, "numpy")
xs = [Link](float(start), float(end), 400)
ys = f(xs)
[Link]()
[Link](xs, ys)
[Link](title)
[Link]()
[Link]()

Purpose: Converts a SymPy expression into a plottable NumPy function and displays its
graph over a given interval.

• [Link](var, expr, "numpy"): This is the bridge between SymPy and NumPy.
SymPy expressions are symbolic objects — they cannot directly process arrays.
lambdify converts the expression into a Python lambda function that uses NumPy
operations, making it compatible with array inputs.
• [Link](start, end, 400): Generates 400 evenly spaced x-values between start
and end. More points = smoother curve. 400 is a good balance between
smoothness and speed.
 • f(xs): Applies the function to every x value in the array at once (vectorized). This is
far more efficient than a Python for loop.

CS Concept: This is an example of the Strategy Pattern — the function accepts any
expression and variable, making it reusable for any mathematical function. The float(start)
call is needed because SymPy symbolic constants like [Link] cannot be directly passed to
NumPy without conversion.

2.3 trap(f, start, end, n) — Trapezoidal Rule

def trap(f, start, end, n):
h = (end - start) / n
s = (f(start) + f(end)) / 2
for i in range(1, n):
s += f(start + i*h)
return h * s

Purpose: Numerically approximates a definite integral using the Trapezoidal Rule.

Mathematical Theory: The Trapezoidal Rule divides the interval [a, b] into n equal
subintervals of width h = (b-a)/n. Instead of computing the exact area under the curve, it
approximates the curve as a series of straight lines connecting adjacent sample points,
forming trapezoids.

T_n = (h/2) * [f(x_0) + 2*f(x_1) + 2*f(x_2) + ... + 2*f(x_{n-

1}) + f(x_n)]

The code equivalent: (f(start) + f(end)) / 2 initializes with the half-weighted endpoints. The
loop adds f(x_i) for each interior point (which corresponds to the coefficient 2 in the formula,
since each interior point is shared by two trapezoids). Finally, multiply by h.

Error Analysis: The error of the Trapezoidal Rule is O(h^2), meaning if you double n
(halving h), the error decreases by a factor of 4. This makes it a second-order method.

CS Concept: f is passed as a function object (first-class function). The parameter f accepts

any callable — a lambda, a def, or a built-in. This is fundamental to functional programming
style in Python.

2.4 simpson(f, start, end, n) — Simpson's Rule

def simpson(f, start, end, n):
if n % 2 != 0:
return "invalid"
h = (end - start) / n
s = f(start) + f(end)
for i in range(1, n):
if i % 2 == 0:
s += 2 * f(start + i*h)
else:
 s += 4 * f(start + i*h)
return h * s / 3

Purpose: Numerically approximates a definite integral using Simpson's Rule, which fits
parabolas through groups of 3 consecutive points instead of straight lines.

Mathematical Theory: Simpson's Rule applies the formula over pairs of subintervals (which
is why n must be even). For each group of 3 points (x_{2k}, x_{2k+1}, x_{2k+2}), it fits an
exact parabola and integrates it analytically. The coefficients in the formula follow the
pattern 1, 4, 2, 4, 2, ..., 4, 1.

S_n = (h/3) * [f(x_0) + 4f(x_1) + 2f(x_2) + 4f(x_3) + ... +

4f(x_{n-1}) + f(x_n)]

In the code: odd-indexed points get coefficient 4, even interior points get coefficient 2, and
the endpoints get coefficient 1. The n % 2 != 0 check enforces the mathematical
requirement that n must be even — an odd n cannot form complete pairs of subintervals.

Error Analysis: Simpson's Rule has error O(h^4), making it a fourth-order method. For the
same n, it is dramatically more accurate than the Trapezoidal Rule. Importantly, Simpson's
Rule gives the EXACT result for polynomials of degree 3 or less, because a parabola can
integrate cubics exactly.

CS Concept: The guard clause if n % 2 != 0: return "invalid" is an early exit pattern. Instead
of wrapping the entire function in an else block, we exit immediately when the precondition
fails. This keeps the main logic unindented and easier to read.

2.5 midpoint_avg(f, start, end, n) — Average Value

def midpoint_avg(f, start, end, n):
h = (end - start) / n
s = 0
for i in range(n):
s += f(start + (i + 0.5) * h)
return s / n

Purpose: Computes the numerical average value of a function over [a, b] by sampling the
function at the midpoint of each subinterval, then dividing by n.

Mathematical Theory: The average value of a continuous function f over [a, b] is:
f_avg = (1 / (b-a)) * integral from a to b of f(x) dx

The code approximates this integral using the Midpoint Rule — for each subinterval, it
samples the function at the center (i + 0.5)*h offset from start. Dividing the sum by n gives
the average. This is equivalent to dividing the integral approximation by (b-a) since h*n =
(b-a).
CS Concept: The expression start + (i + 0.5) * h cleverly computes the midpoint of the i-th
subinterval without needing a separate variable. This is a common trick in numerical
computing to avoid off-by-one errors.
3. Exercise-by-Exercise Explanation
Exercise 1 — Fundamental Definite Integrals
e1 = [("a", [Link](x**3 + 2*x**2 + 3, (x, 1, 2))), ...]
show("Exercise 1", e1)

This exercise is a collection of 17 definite integrals, each testing a different integration

technique. They all use the Fundamental Theorem of Calculus:
integral from a to b of f(x) dx = F(b) - F(a), where F'(x) =
f(x)

SymPy's [Link](expr, (x, a, b)) applies this theorem symbolically — it first finds the
antiderivative F(x), then evaluates F(b) - F(a) exactly. The results are rational numbers or
expressions involving pi, log, sqrt, etc.

Key techniques practiced across the 17 integrals:

# Technique Example / Notes

a–f Power rule integral of x^n dx = x^(n+1)/(n+1) + C. Works for
fractional powers too (e.g., x*sqrt(x) = x^(3/2)).
g, Trig identities Uses identities like 1/sin^2(2x) = csc^2(2x), and the
h identity sin(2x)*cos(2x) = sin(4x)/2 to simplify before
integrating.
i, j Exponential integrals Expands products involving e^x, then integrates each
term. Note: integral of e^x dx = e^x + C.
k Base-a exponential integral of 2^x dx = 2^x / ln(2) + C. SymPy handles this
automatically.
m Partial fractions 1/(x*(x+1)) = 1/x - 1/(x+1). Each term is then integrated
as a logarithm.
n, Absolute value integrals SymPy splits the interval at the zero of the expression
o inside abs(), integrating each piece separately.
p Square root (completing sqrt(x^2 - 3x + 2) = sqrt((x-1)(x-2)) requires interval
square) splitting since the expression is zero at x=1 and x=2.
q Trig simplification sqrt(1 + cos(2x)) = sqrt(2) * |cos(x)|. The result over [0,
pi] involves careful handling of the sign change.

Why use SymPy and not just evaluate numerically? SymPy gives exact answers like log(2),
sqrt(2)*pi/4, or 13/30 — not floating-point approximations like 0.6931... This is crucial when
precision matters, such as in proofs or further symbolic computation.
Exercise 2 — Advanced Integrals and Double Integral
e2 = [
("a", [Link](x**3 - 3*[Link](x)*[Link](x), (x, 0, [Link]/2))),
("b", [Link]([Link](x**2)**2, (x, 0, 1))),
("c", [Link]([Link](1 + (x**2+1) + (x+1)**2), (x, 0, 3))),
("d", [Link]([Link](x**2*y, (x, 0, 3)), (y, 1, 2)))
]

(a) The integrand x^3 - 3*sin(x)*cos(x) is split linearly. The key identity is: sin(x)*cos(x) =
sin(2x)/2. So the second term becomes -3/2 * integral of sin(2x) dx, which equals 3/4 *
cos(2x) + C.

(b) sin^2(x^2) has no elementary antiderivative — it cannot be expressed in terms of basic

functions. SymPy returns a result in terms of special functions (Fresnel integrals). This is a
good example of why symbolic math tools are needed: a calculator would only give a
number, not the form.

(c) The expression under the square root simplifies to sqrt(x^4 + 4x^2 + 2x + 3). This form
appears in arc-length computations: the arc length of a curve y=f(x) from a to b is integral of
sqrt(1 + [f'(x)]^2) dx.

(d) Double Integral: This uses Fubini's Theorem — for a rectangular region, a double
integral can be evaluated as two successive single integrals:
integral from 1 to 2 (integral from 0 to 3 of x^2*y dx) dy

The code nests [Link]() calls, exactly mirroring the notation. The inner integral treats
y as a constant: integral of x^2*y dx = y * x^3/3. Evaluated at x=0 and x=3 gives 9y. Then
the outer integral: integral from 1 to 2 of 9y dy = 9 * [y^2/2] = 9*(4/2 - 1/2) = 13.5.

The three plot_expr() calls after this exercise display graphs of (a), (b), and (c) — helping
visualize what you just integrated.

CS Concept: The nested [Link]() calls demonstrate function composition. The result of
the inner call (an expression in y) is passed directly as the expression argument to the outer
call. This mirrors mathematical notation precisely and is an example of clean, readable
code.

Exercise 3 — Average Value of a Function

for k, expr, start, end in e3:
avg = [Link](expr, (x, start, end)) / (end - start)
print(f"({k}) {[Link](avg)}")
plot_expr(expr, x, start, end, "Exercise 3" + k)

This exercise applies the Mean Value Theorem for Integrals. The average value of a
continuous function f over [a, b] is:
f_avg = (1 / (b - a)) * integral from a to b of f(x) dx
Geometrically, f_avg is the height of a rectangle with base (b-a) whose area equals the
area under the curve. The code computes the integral symbolically with SymPy, then
divides by (end - start). This gives an exact symbolic result, which is then simplified.

The four sub-problems cover:

• (a) f(x) = x^2 - 1 on [0, sqrt(3)]: Average of a parabola. The result is a rational
number.
• (b) f(x) = -x^2/2 on [0, 3]: Downward parabola. The average is negative since the
function is non-positive.
• (c) f(x) = -3x^2 - 1 on [0, 1]: Always negative — average confirms this.
• (d) f(x) = x^2 - x on [-2, 1]: Mixed sign function. The average captures the net
balance of positive and negative regions.

Each sub-problem also plots the function, so you can visually verify that the average value
looks correct relative to the graph.

CS Concept: The for loop unpacks a 4-element tuple (k, expr, start, end) from the list e3.
This is a clean pattern for iterating over structured data. Each tuple stores all the information
needed for one sub-problem, keeping the code organized and avoiding repeated variable
names.

Exercise 4 — Improper Integral

e4a = [Link](x**2 * [Link](x), (x, -4, 9))
e4b = [Link]([Link](-[Link](1, 2) * x**2), (x, -[Link], [Link]))

(a) A standard definite integral using integration by parts twice (since x^2 differentiates to 0
after two steps, and cos(x) cycles). SymPy handles this automatically. The result is an
exact expression involving sin and cos evaluated at -4 and 9.

(b) The Gaussian Integral: This is one of the most famous improper integrals in all of
mathematics:
integral from -infinity to +infinity of e^(-x^2/2) dx =
sqrt(2*pi)

It cannot be evaluated by finding an elementary antiderivative — there is none. Instead, the

proof uses a clever polar coordinates trick (squaring the integral and converting to 2D).
SymPy uses its knowledge of special functions to return the exact value sqrt(2*pi).

This integral is the foundation of the normal distribution in probability and statistics. The
standard normal PDF is (1/sqrt(2*pi)) * e^(-x^2/2), which integrates to exactly 1 over all
reals.

[Link](1, 2) creates the exact fraction 1/2 as a SymPy object. Using 0.5 (a float) would
introduce floating-point representation into the symbolic calculation, potentially giving less
clean results. Always prefer [Link] for exact fractions in symbolic math.
 [Link] is SymPy's representation of infinity. Passing -[Link] and [Link] as the limits of
integration tells SymPy this is an improper integral and triggers the appropriate limit
evaluation.

Exercise 5 — Physics Application: Velocity to Displacement

v = 160 - 32*t
print([Link](v, (t, 0, 8)))

This models a ball thrown straight up with initial velocity 160 ft/s under Earth's gravity (g =
32 ft/s^2). The velocity function is v(t) = 160 - 32t.

Key physics: displacement equals the integral of velocity over time:

displacement = integral from 0 to 8 of (160 - 32t) dt

Evaluating: [160t - 16t^2] from 0 to 8 = (160*8 - 16*64) - 0 = 1280 - 1024 = 256 feet.

Note that this is NET displacement (final position minus initial position), not total distance
traveled. The ball reaches its peak when v(t) = 0, i.e., at t = 5 seconds. After that, it falls
back down. The integral from 0 to 8 gives the net height above the starting point at t=8.

CS/Physics Connection: v = 160 - 32*t defines t as a SymPy symbol (declared at the top: t =
[Link]('t')). The expression is symbolic — Python doesn't compute a number when you
write this, it creates a symbolic expression object. Integration then happens when
[Link]() is called.

Exercise 6 — Economics Application: Marginal Cost

dc = 1 / (2 * [Link](x))
print([Link](dc, (x, 1, 100)))

In economics, if C(x) is the total cost of producing x units, then C'(x) (the derivative) is the
marginal cost — the cost of producing one more unit. Here, the marginal cost function is:
C'(x) = 1 / (2 * sqrt(x))

Integrating the marginal cost from x=1 to x=100 gives the change in total cost as production
increases from 1 to 100 units:
Delta C = integral from 1 to 100 of 1/(2*sqrt(x)) dx =
[sqrt(x)] from 1 to 100 = 10 - 1 = 9

The antiderivative of 1/(2*sqrt(x)) = (1/2)*x^(-1/2) is x^(1/2) = sqrt(x). So the total additional

cost is exactly 9 units (dollars, in whatever currency the problem uses).

This is a perfect example of the Fundamental Theorem of Calculus in an applied context:

the definite integral of a rate of change (marginal cost) gives the total change (cost
increase).
Exercise 7 — Growth Modeling Over Time
H = [Link](t + 1) + 5*t**[Link](1, 3)
print("(a)", [Link](t, 0), [Link](t, 4), [Link](t, 8))
print("(b)", [Link]([Link](H, (t, 0, 8)) / 8))

H(t) = sqrt(t+1) + 5*t^(1/3) models a quantity that grows over time — for example, the
height of a plant in centimeters after t hours.

Part (a) — Point evaluations: [Link](t, 0) substitutes t=0 into the expression symbolically.
subs() is SymPy's substitution method. The three values at t=0, 4, 8 give snapshots of
growth at specific times.

Part (b) — Time-average: Computes the average value of H over 8 hours using the Mean
Value Theorem for Integrals:
H_avg = (1/8) * integral from 0 to 8 of H(t) dt

[Link](1, 3) is used instead of the float 1/3 to ensure exact computation of t^(1/3). In
Python, 1/3 = 0.3333... (a float), which introduces rounding. [Link](1, 3) is the exact
fraction 1/3 in symbolic form.

Exercise 8 — Numerical Average Value with Convergence Testing

e8 = [
("a", lambda z: 1 - z, 0, 1),
("b", lambda z: z**2 + 1, 0, 1),
("c", lambda z: [Link](z), -[Link], [Link]),
("d", lambda z: abs(z), -1, 1)
]
for k, f, start, end in e8:
for n in [4, 100, 200, 1000]:
print(n, midpoint_avg(f, start, end, n))

This exercise tests how the numerical average converges as n increases. Lambda
functions are used here instead of def — they are anonymous one-line functions
appropriate for simple expressions.

The nested for loop runs each function at n = 4, 100, 200, 1000 subdivisions. You should
observe that the values stabilize (converge) as n grows, demonstrating the fundamental
principle: more subdivisions = more accurate approximation.

Expected convergence behavior:

• (a) 1 - z on [0, 1]: Exact average = 0.5. Even at n=4, the midpoint rule is exact for
linear functions.
• (b) z^2 + 1 on [0, 1]: Exact average = 4/3 ≈ 1.333... Converges quickly.
 • (c) cos(z) on [-pi, pi]: Exact average = 0 by symmetry. Note [Link] is used (not
[Link]) since these are NumPy lambda functions.
• (d) |z| on [-1, 1]: Exact average = 0.5. The non-smooth point at z=0 slows
convergence slightly.

CS Concept: lambda z: 1 - z defines an anonymous function inline. It is equivalent to def

f(z): return 1 - z but shorter. Lambdas are commonly used in Python when passing simple
functions as arguments — a core functional programming technique.

Exercise 9 — Trapezoidal Rule: Convergence Study

print(trap(lambda z: [Link](-z**2), 0, 1, 3))
for n in [1, 3, 4, 6]:
print(n, trap(lambda z: 2*z**2 + 5*z + 12, -1, 5, n))

Exercise 9 applies the trap() function to four different integrands at increasing values of n,
demonstrating convergence and testing the accuracy limits of the trapezoidal method.

(a) e^(-z^2) on [0, 1]: This function has no elementary antiderivative — it is the basis of the
error function erf(x). Numerical methods are the only practical way to evaluate it. At n=3,
the result is an approximation; increasing n improves accuracy.

(b) 2z^2 + 5z + 12 on [-1, 5] at n=1,3,4,6: A degree-2 polynomial. The trapezoidal rule is

NOT exact for quadratics (it is only exact for linear functions). You will see the result
converging toward the true value as n increases. The exact integral can be computed by
hand: [(2/3)z^3 + (5/2)z^2 + 12z] from -1 to 5 = a specific rational number.

(c) z^3 + 2z^2 - 5z - 2 on [0, 2]: A cubic polynomial. Tested at n=2,4,6,8. The trapezoidal
rule converges, but more slowly than Simpson's Rule for the same function.

(d) z*e^(-z) on [0.2, 3.8]: Requires integration by parts to solve analytically: integral of
z*e^(-z) dz = -e^(-z)(z+1) + C. The trapezoidal approximation converges well since the
function is smooth.

CS Concept: The for n in [1, 3, 4, 6] pattern iterates over specific values (not a range). This
is intentional — the exercise checks specific n values that may be required by a textbook
problem, not a smooth sequence. Python lists can hold any arbitrary values for this purpose.

Exercise 10 — Simpson's Rule: Convergence Study

print(simpson(lambda z: [Link](-z**2), 0, 1, 3))
for n in [1, 3, 4, 6]:
print(n, simpson(lambda z: 2*z**2 + 5*z + 12, -1, 5, n))

Exercise 10 repeats the same four integrands from Exercise 9, but uses the simpson()
function instead. This allows direct comparison of accuracy between the two methods.
Key comparison points:
• (a) e^(-z^2) at n=3: simpson() returns 'invalid' because 3 is odd (Simpson's Rule
requires even n). This is caught by the n % 2 != 0 guard clause.
• (b) Quadratic at n=4,6: Simpson's Rule gives the EXACT answer for quadratic (and
cubic) polynomials, because a parabola integrates degree-3 polynomials exactly.
You will see the result match the true value precisely, even at small n.
• (c,d) Higher n: For n=4,6,8, Simpson's Rule converges significantly faster than the
Trapezoidal Rule. The difference in accuracy at the same n shows the practical
advantage of Simpson's fourth-order method.

Property Trapezoidal Rule Midpoint Rule Simpson's Rule

Shape used Trapezoids (linear) Rectangles Parabolas
(midpoint) (quadratic)
Coefficients 1, 2, 2, ..., 2, 1 1, 1, 1, ..., 1 1, 4, 2, 4, ..., 4, 1
Error order O(h^2) O(h^2) O(h^4)
n restriction Any n >= 1 Any n >= 1 n must be even
Exact for Linear polynomials Linear polynomials Polynomials degree
<= 3
4. Key CS and Programming Patterns
Beyond the mathematics, this code demonstrates several important programming concepts
that are valuable for a Computer Science student's portfolio and career.

Pattern / Concept Where Used Why It Matters

First-class trap(f, ...), simpson(f, ...) — f Core to functional programming;
functions is passed as an argument allows reusable, flexible algorithms
Lambda functions lambda z: [Link](-z**2) in Concise inline functions; commonly
Exercise 8–10 used with map(), filter(), sorted()
Tuple unpacking for k, expr, start, end in e3: Clean pattern for iterating structured
data; avoids indexing
Guard clauses if n % 2 != 0: return 'invalid' Early exit pattern; keeps main logic
clean and readable
Library [Link](), [Link]() Using black-box APIs; understanding
abstraction input/output contracts
Vectorized ys = f(xs) where xs is a Processing entire arrays at once;
operations NumPy array core to scientific Python (NumPy)
List of tuples e1 = [("a", expr_a), ("b", Organizing related data; pattern
(data structure) expr_b), ...] appears in databases, CSV parsing
Convergence for n in [4, 100, 200, 1000]: Empirical validation of algorithms;
testing print(midpoint_avg(...)) standard in numerical computing
5. Library Quick Reference
5.1 SymPy Cheatsheet
Function / Object Description
[Link]('x y t a') Declare symbolic variables — must be done before use
[Link](expr, x) Indefinite integral (antiderivative) of expr with respect to x
[Link](expr, (x, Definite integral from a to b
a, b))
[Link](expr) Simplify expression to a shorter/cleaner form
[Link](x, value) Substitute a value for a symbol in an expression
[Link](x, expr, Convert SymPy expression to a NumPy-compatible
'numpy') function
[Link](p, q) Exact fraction p/q (avoids floating-point errors)
[Link] Infinity (used as integration limit for improper integrals)
[Link], [Link](x), Exact constants and functions
[Link](x)

5.2 NumPy Cheatsheet

Function Description
[Link](a, b, n) Create n evenly spaced values between a and b
(inclusive)
[Link](x) Element-wise e^x for arrays
[Link](x), [Link](x) Element-wise trig functions for arrays
float(sym_value) Convert SymPy result to Python float for NumPy
compatibility

5.3 Matplotlib Cheatsheet

Function Description
[Link]() Create a new plot figure
[Link](xs, ys) Plot y-values against x-values as a line graph
[Link](text) Set the plot title
Function Description
[Link]() Add a background grid for readability
[Link]() Render and display the plot
6. How to Run This Code
6.1 Setup
Install the required libraries using pip:
pip install sympy numpy matplotlib

Then run the script from the terminal:

python calculus_toolkit.py

6.2 Expected Output

The terminal will print all numerical results for Exercises 1–10. Matplotlib will open separate
plot windows for each exercise that calls plot_expr(). Close each window to continue to the
next plot.

6.3 Tips for Experimentation

• Change the limits: Modify (x, 1, 2) to different values and observe how results
change.
• Increase n in numerical methods: Try n=10000 in trap() or simpson() and see how
much the result changes.
• Try a new integral: Add a new entry to e1 with any expression you want to
integrate.
• Compare symbolic vs numerical: For any function, compare [Link]() (exact)
with trap() or simpson() (approximate) to measure the error.

Generated as a study reference — Python Calculus Toolkit Explanation