Calculus with Python: SymPy and NumPy Guide

What Is Calculus With Python?

Calculus with Python is the practice of computing limits, derivatives, integrals, series and differential equations with Python libraries instead of (or alongside) pencil and paper. The work splits into two kinds. Symbolic calculus, done with SymPy, returns exact formulas such as $\frac{x^3}{3}$ or $\sqrt{\pi}$. Numerical calculus, done with NumPy (and SciPy where it is installed), returns decimal approximations such as $0.7468$.

The distinction matters because the two answer different questions. A symbolic result tells you the rule; a numerical result tells you a value. Calculus with Python works best when you use each where it fits:

Every block of code below was run with Python 3.14, SymPy 1.14 and NumPy 2.5, and the output shown is the real output.

How Do You Do Calculus With Python?

Every calculus task in SymPy follows the same four moves:

1. Import the functions you need from sympy.

2. Declare symbols with symbols('x'), so Python knows x is a variable, not a number.

3. Build the expression with ** for powers and SymPy's own sin, exp and log.

4. Call the operation (limit, diff, integrate, series, dsolve) and print the result.

Step 2 is the one beginners skip. Without it, x is an undefined name and Python stops with a NameError before any calculus happens. Installing the libraries is a single pip install sympy numpy in a terminal; add scipy and matplotlib when you need them.

How Do You Compute Limits In SymPy?

A limit in SymPy is computed with limit(f, x, a), which returns the exact value $f(x)$ approaches as $x \to a$. Infinity is written oo (two lowercase o's), and the optional dir argument picks a one-sided limit. This is the same idea as the limit of a function on paper:

$$\lim_{x \to 0} \frac{\sin x}{x} = 1, \qquad \lim_{x \to \infty} \left(1 + \frac{1}{x}\right)^{x} = e$$

from sympy import symbols, limit, sin, oo

x = symbols('x')
print(limit(sin(x)/x, x, 0))
print(limit((1 + 1/x)**x, x, oo))
print(limit(1/x, x, 0, dir='+'))
print(limit(1/x, x, 0, dir='-'))
print(limit(1/x, x, 0, dir='+-'))

1
E
oo
-oo
zoo

Read the last three lines carefully. From the right, $\frac{1}{x}$ climbs to $+\infty$; from the left it falls to $-\infty$. When you ask for both sides with dir='+-', SymPy answers zoo (complex infinity), its way of saying the two-sided limit does not exist. That name looks like a typo the first time you meet it.

SymPy's default direction is '+', so a bare limit(1/x, x, 0) quietly reports only the right-hand limit.

How Do You Find Derivatives In Python?

A derivative in SymPy is one call: diff(f, x) for the first derivative and diff(f, x, n) for the $n$th. The result is a new SymPy expression, so you can simplify it, substitute into it or differentiate it again. This is the derivative of a function computed by the same rules you learn by hand (product, chain, quotient), applied without arithmetic errors.

from sympy import symbols, diff, sin, exp, pi

x = symbols('x')
f = x**3 * sin(x)
print(diff(f, x))
print(diff(f, x, 2))
print(diff(exp(x**2), x))
print(diff(f, x).subs(x, pi))

x**3*cos(x) + 3*x**2*sin(x)
x*(-x**2*sin(x) + 6*x*cos(x) + 6*sin(x))
2*x*exp(x**2)
-pi**3

The first line is the product rule, $\frac{d}{dx}\left[x^3 \sin x\right] = x^3 \cos x + 3x^2 \sin x$. The third is the chain rule: the inner derivative $2x$ multiplies $e^{x^2}$. The last line evaluates $f'(\pi) = \pi^3 \cos \pi + 3\pi^2 \sin \pi = -\pi^3$ exactly, with no decimal in sight.

"Is it possible to do abstract differentiation and return a function instead of just a numerical value?" That question, asked on a Physics Forums thread about a home-made derivative calculator, has a direct answer: yes, and SymPy is the tool. diff returns a formula. When you then need fast numbers from that formula, lambdify turns it into an ordinary Python function:

from sympy import symbols, diff, sin, lambdify

x = symbols('x')
df = diff(x**3 * sin(x), x)
df_num = lambdify(x, df, 'numpy')
print(df_num(1.0))

3.064715260291829

That value is $\cos 1 + 3\sin 1 \approx 3.0647$, the slope of $x^3 \sin x$ at $x = 1$.

For functions of several variables, pass the variable you differentiate with respect to. Every other symbol is held constant, which is exactly a partial derivative:

from sympy import symbols, diff

x, y = symbols('x y')
g = x**2 * y + 3*x*y**3
print(diff(g, x))
print(diff(g, y))
print(diff(g, x, y))

2*x*y + 3*y**3
x**2 + 9*x*y**2
2*x + 9*y**2

The mixed partial $\frac{\partial^2 g}{\partial x,\partial y} = 2x + 9y^2$ comes out the same whichever order you differentiate in, as Clairaut's theorem promises for smooth functions.

How Do You Compute A Derivative Numerically?

A numerical derivative approximates the slope from function values alone, using a difference quotient with a small step $h$. The central version is the most accurate simple choice:

$$f'(a) \approx \frac{f(a + h) - f(a - h)}{2h}$$

Here $f$ is any function that takes a number and returns a number, $a$ is the point, and $h$ is a small positive step.

import numpy as np

def derivative(f, a, h=1e-5):
    return (f(a + h) - f(a - h)) / (2 * h)

print(derivative(np.sin, 0.0))
print(derivative(np.exp, 1.0))

0.9999999999833332
2.718281828517632

The exact answers are $\cos 0 = 1$ and $e^1 \approx 2.718281828$, so both are correct to about ten digits. A smaller $h$ is not automatically better, though. Watch the forward difference $\frac{\sin(1 + h) - \sin 1}{h}$ against the exact $\cos 1$:

import numpy as np

exact = np.cos(1.0)
for h in [1e-1, 1e-3, 1e-5, 1e-8, 1e-11, 1e-14]:
    approx = (np.sin(1.0 + h) - np.sin(1.0)) / h
    print(f"h = {h:.0e}   error = {abs(approx - exact):.2e}")

h = 1e-01   error = 4.29e-02
h = 1e-03   error = 4.21e-04
h = 1e-05   error = 4.21e-06
h = 1e-08   error = 1.41e-08
h = 1e-11   error = 1.17e-06
h = 1e-14   error = 3.71e-03

The error shrinks until about $h = 10^{-8}$, then grows again. Below that point the subtraction $\sin(1 + h) - \sin 1$ cancels almost every stored digit, and floating-point round-off takes over. This is the numerical version of the limit definition, and it shows why the limit is a mathematical idea that a computer can only approach.

How Do You Find Integrals In Python?

SymPy's integrate handles both kinds of integral. integrate(f, x) returns an antiderivative, and integrate(f, (x, a, b)) returns the exact value of a definite integral. Infinite bounds use oo, which covers improper integrals too.

from sympy import symbols, integrate, exp, oo

x = symbols('x')
print(integrate(x**2, x))
print(integrate(x * exp(x), x))
print(integrate(x**2, (x, 0, 3)))
print(integrate(exp(-x), (x, 0, oo)))

x**3/3
(x - 1)*exp(x)
9
1

In standard notation those four lines say:

$$\int x^2,dx = \frac{x^3}{3} + C, \qquad \int x e^{x},dx = (x - 1)e^{x} + C$$

$$\int_{0}^{3} x^2,dx = 9, \qquad \int_{0}^{\infty} e^{-x},dx = 1$$

Notice what SymPy left out. It prints x**3/3, not x**3/3 + C, because it returns one antiderivative, not the whole family. You add the constant yourself whenever you write the answer down. Differentiating back confirms each result: $\frac{d}{dx}\left[(x - 1)e^{x}\right] = e^{x} + (x - 1)e^{x} = x e^{x}$.

Some integrals have no formula built from polynomials, roots, exponentials, logarithms and trig functions. SymPy says so plainly, by answering with a named special function:

from sympy import symbols, integrate, exp, sin

x = symbols('x')
print(integrate(exp(-x**2), x))
area = integrate(sin(x)/x, (x, 0, 1))
print(area)
print(area.evalf())

sqrt(pi)*erf(x)/2
Si(1)
0.946083070367183

$\int e^{-x^2},dx$ is non-elementary, so SymPy writes it with the error function, $\frac{\sqrt{\pi}}{2}\operatorname{erf}(x) + C$, defined in the NIST Digital Library of Mathematical Functions. Likewise $\int_{0}^{1} \frac{\sin x}{x},dx$ comes back as the sine integral $\operatorname{Si}(1)$, and .evalf() turns it into $0.9461$ (to 4 decimal places).

How Do You Integrate Numerically With NumPy?

When a function has no usable antiderivative, or you only have data points, you integrate numerically. NumPy's np.trapezoid (new in NumPy 2.0, replacing the older np.trapz) applies the trapezoidal rule to sampled values:

import numpy as np

x = np.linspace(0, 1, 101)
y = np.exp(-x**2)
print(np.trapezoid(y, x))

0.7468180014679698

The exact value is $\int_{0}^{1} e^{-x^2},dx = \frac{\sqrt{\pi}}{2}\operatorname{erf}(1) \approx 0.746824132812427$ (SymPy, 15 digits). One hundred trapezoids land within about $6 \times 10^{-6}$. Fitting parabolas instead of straight lines, the idea behind Simpson's rule, does far better with the same 101 points:

import numpy as np

def simpson(f, a, b, n=100):
    x = np.linspace(a, b, n + 1)
    y = f(x)
    h = (b - a) / n
    return h / 3 * (y[0] + 4 * y[1:-1:2].sum() + 2 * y[2:-1:2].sum() + y[-1])

print(simpson(lambda t: np.exp(-t**2), 0, 1))

0.7468241328941763

That matches the exact value to about $8 \times 10^{-11}$. If SciPy is installed, scipy.integrate.quad(f, a, b) is the standard adaptive routine; per the SciPy documentation it returns two numbers, the integral and an estimate of its absolute error. The ideas behind all three methods are covered in numerical integration.

How Do You Get Taylor Series In Python?

series(f, x, a, n) expands $f$ around $x = a$ up to (but not including) the power $x^n$, and appends an order term O(x**n) that marks what was dropped. removeO() strips that term so you can evaluate the polynomial. These are the Taylor series you derive by hand from repeated derivatives.

from sympy import symbols, series, exp, sin

x = symbols('x')
print(series(exp(x), x, 0, 6))
print(series(sin(x), x, 0, 8).removeO())

1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + O(x**6)
-x**7/5040 + x**5/120 - x**3/6 + x

The denominators $2, 6, 24, 120, 5040$ are the factorials of $2, 3, 4, 5$ and $7$, the pattern in

$$e^{x} = 1 + x + \frac{x^2}{2} + \frac{x^3}{6} + \frac{x^4}{24} + \frac{x^5}{120} + \cdots$$

SymPy prints the sine polynomial with the highest power first; the order of terms in the output is SymPy's display choice, not a different answer.

How Do You Solve Differential Equations In Python?

dsolve solves ordinary differential equations. You declare the unknown with Function('y'), write the equation with Eq, and optionally pass initial conditions through ics. Without conditions you get the general solution with constants C1, C2; with them you get the particular solution of an initial value problem.

from sympy import symbols, Function, Eq, dsolve

t = symbols('t')
y = Function('y')
ode = Eq(y(t).diff(t), -2 * y(t))
print(dsolve(ode, y(t)))
print(dsolve(ode, y(t), ics={y(0): 5}))

Eq(y(t), C1*exp(-2*t))
Eq(y(t), 5*exp(-2*t))

For $\frac{dy}{dt} = -2y$ the general solution is $y = C_1 e^{-2t}$, and $y(0) = 5$ fixes $C_1 = 5$. Second-order equations work the same way, and checkodesol substitutes the answer back for you:

from sympy import symbols, Function, Eq, dsolve, checkodesol

t = symbols('t')
y = Function('y')
ode = Eq(y(t).diff(t, 2) + y(t), 0)
sol = dsolve(ode, y(t))
print(sol)
print(checkodesol(ode, sol))

Eq(y(t), C1*sin(t) + C2*cos(t))
(True, 0)

(True, 0) means the residual after substitution is zero: $y'' + y = 0$ is satisfied by $C_1 \sin t + C_2 \cos t$ for every $C_1$ and $C_2$.

How Do You Plot Calculus In Python?

Plotting turns a formula into a picture you can check by eye. With Matplotlib, the usual pattern is to sample with np.linspace, draw the curve with plt.plot, and shade the area under it with plt.fill_between(x, y), which per the Matplotlib documentation fills the region between two curves (the second defaults to $y = 0$). For a quick look with no extra install, SymPy's textplot draws the graph as text:

from sympy import symbols
from sympy.plotting import textplot

x = symbols('x')
textplot(x**2 - 2, -2, 2, W=40, H=12)

      2 |\                                      /|
        | \                                    / |
        |  \                                  /  |
        |   \                                /   |
        |    \                              /    |
0.18181 |-----\----------------------------/-----|
        |      ..                        ..      |
        |        \                      /        |
        |         ..                  ..         |
        |           ..              ..           |
        |             ...        ...             |
     -2 |________________________________________|
         -2                  0.0512821          2|

The parabola $y = x^2 - 2$ dips to $-2$ at $x = 0$ and rises to $2$ at $x = \pm 2$. The middle axis labels are textplot's own sampled values, not features of the function.

When Should You Use SymPy And When NumPy?

Use SymPy for exactness and NumPy for speed. SymPy carries exact values (fractions, $\pi$, $\sqrt{2}$) through every step, so integrate(x**2, (x, 0, 3)) returns the integer $9$, not $8.999999$. NumPy works in 64-bit floating point, which stores roughly 15 to 16 significant digits and is far faster on large arrays.

The practical rule:

• Use SymPy when you need a formula, an exact value, or a check on hand work.

• Use NumPy (or SciPy) when you need numbers fast, from data or from an integral with no closed form.

• Use both by deriving symbolically, then calling lambdify to evaluate numerically.

Examples Of Calculus With Python

Example 1

Find the derivative of $f(x) = x^2 e^{x}$ and evaluate it at $x = 1$.

from sympy import symbols, diff, exp, factor

x = symbols('x')
f = x**2 * exp(x)
df = diff(f, x)
print(factor(df))
print(df.subs(x, 1))
print(df.subs(x, 1).evalf(6))

x*(x + 2)*exp(x)
3*E
8.15485

$$f'(x) = 2x e^{x} + x^2 e^{x}$$

$$f'(x) = x(x + 2)e^{x}$$

$$f'(1) = 1 \cdot 3 \cdot e = 3e$$

Final answer: $f'(x) = x(x + 2)e^{x}$ and $f'(1) = 3e \approx 8.1548$.

Example 2

Estimate $\lim_{n \to \infty} \left(1 + \frac{1}{n}\right)^{n}$ with Python.

The tempting route is to plug in bigger and bigger numbers with ordinary floats:

for n in [10, 1000, 10**8, 10**16]:
    print(n, (1 + 1/n)**n)

10 2.5937424601000023
1000 2.7169239322355936
100000000 2.7182817983473577
10000000000000000 1.0

The values creep toward $2.71828$ and then collapse to exactly $1.0$. That cannot be right: every term $\left(1 + \frac{1}{n}\right)^{n}$ is bigger than $2$ for $n \geq 2$. The break is round-off. For $n = 10^{16}$, the number $\frac{1}{n}$ is smaller than the gap between $1$ and the next representable float, so 1 + 1/n is stored as exactly 1.0, and $1^{n} = 1$.

The rescue is to ask for the limit symbolically, so no float ever appears:

from sympy import symbols, limit, oo

n = symbols('n', positive=True)
L = limit((1 + 1/n)**n, n, oo)
print(L, L.evalf(15))

E 2.71828182845905

Final answer: the limit is $e \approx 2.7183$. Numerical plugging gives evidence; limit gives the answer.

Example 3

Find the area enclosed between $y = x$ and $y = x^2$.

The curves meet where $x = x^2$. On $[0, 1]$ the line sits above the parabola, so the area is the integral of top minus bottom, exactly as in area between two curves.

from sympy import symbols, solve, integrate

x = symbols('x')
print(solve(x - x**2, x))
print(integrate(x - x**2, (x, 0, 1)))

[0, 1]
1/6

$$A = \int_{0}^{1} \left(x - x^2\right),dx$$

$$A = \left[\frac{x^2}{2} - \frac{x^3}{3}\right]_{0}^{1}$$

$$A = \frac{1}{2} - \frac{1}{3} = \frac{1}{6}$$

Final answer: $\frac{1}{6}$ square units.

Example 4

Find and classify the critical points of $f(x) = x^3 - 6x^2 + 9x + 1$.

from sympy import symbols, diff, solve

x = symbols('x')
f = x**3 - 6*x**2 + 9*x + 1
crit = solve(diff(f, x), x)
for c in crit:
    print(c, diff(f, x, 2).subs(x, c), f.subs(x, c))

1 -6 5
3 6 1

Each line prints the critical point, $f''$ there, and $f$ there.

$$f'(x) = 3x^2 - 12x + 9 = 3(x - 1)(x - 3)$$

$$f''(1) = -6 < 0, \qquad f''(3) = 6 > 0$$

Final answer: a local maximum at $(1, 5)$ and a local minimum at $(3, 1)$, the same verdict the second derivative test gives for maxima and minima by hand.

Example 5

Approximate $\sin\left(\frac{1}{2}\right)$ with the degree-5 Taylor polynomial and measure the error.

from sympy import symbols, series, sin, Rational

x = symbols('x')
p = series(sin(x), x, 0, 6).removeO()
approx = p.subs(x, Rational(1, 2))
print(p)
print(approx, approx.evalf(10))
print(sin(Rational(1, 2)).evalf(10))
print((approx - sin(Rational(1, 2))).evalf(3))

x**5/120 - x**3/6 + x
1841/3840 0.4794270833
0.4794255386
1.54e-6

$$p(x) = x - \frac{x^3}{6} + \frac{x^5}{120}$$

$$p\left(\tfrac{1}{2}\right) = \frac{1}{2} - \frac{1}{48} + \frac{1}{3840} = \frac{1841}{3840}$$

Final answer: $p\left(\frac{1}{2}\right) \approx 0.4794$, within about $1.54 \times 10^{-6}$ of $\sin\left(\frac{1}{2}\right)$. Using Rational(1, 2) instead of 0.5 keeps the polynomial value exact.

Example 6

A cup of tea at $90^\circ\text{C}$ cools in a $20^\circ\text{C}$ room by Newton's law of cooling, $\frac{dT}{dt} = -0.1,(T - 20)$, with $t$ in minutes. Find $T(t)$ and the temperature after 10 minutes.

from sympy import symbols, Function, Eq, dsolve, Rational

t = symbols('t')
T = Function('T')
k = Rational(1, 10)
ode = Eq(T(t).diff(t), -k * (T(t) - 20))
sol = dsolve(ode, T(t), ics={T(0): 90})
print(sol)
print(sol.rhs.subs(t, 10).evalf(6))

Eq(T(t), 20 + 70*exp(-t/10))
45.7516

$$T(t) = 20 + 70e^{-t/10}$$

$$T(10) = 20 + \frac{70}{e} \approx 45.7516$$

Final answer: $T(t) = 20 + 70e^{-t/10}$, and the tea is about $45.75^\circ\text{C}$ after 10 minutes. The same model, with a different sign, describes exponential growth and decay.

Why Does Calculus With Python Matter Beyond Homework?

Many of the integrals that turn up in real problems have no elementary formula at all, and that single fact is why calculus moved onto computers. Hand methods cover the integrals and equations that were chosen because they can be solved by hand; the ones that turn up in physics, engineering and data work usually cannot.

• Problems with no closed form. $\int_{0}^{1} e^{-x^2},dx$ underlies the normal distribution, yet it has no elementary antiderivative. Numerical integration is the only route to its value, $0.7468$.

• Checking work you did by hand. A derivative typed into diff either matches your page or it does not. That instant feedback is how many learners finally trust (or catch) their product and chain rules.

• Scale. Training a machine-learning model means computing derivatives of a loss function millions of times, which no person does by hand. The gradient in gradient descent is the multivariable cousin of the diff(g, x) call above.

The destination is worth naming: the limit, derivative and integral you learn on paper are the same objects that simulation and machine-learning code compute every second. Python only does the arithmetic; deciding which quantity to compute is still calculus.

Where Do Students Trip Up On Calculus With Python?

Most errors in calculus code are not calculus errors. They are places where Python's rules differ from the notation on paper, and the code runs anyway (or fails with a message that does not mention calculus).

Mistake 1: Writing Powers With ^

Where it slips in: Typing $x^2$ the way a calculator or textbook shows it.

Don't do this:x^2. In Python ^ is bitwise XOR, so 3^2 evaluates to 1, and with a SymPy symbol it fails with TypeError: unsupported operand type(s) for ^: 'Symbol' and 'int'.

The correct way:x**2. Python's power operator is **, as the Python language reference states.

The first instinct of nearly every newcomer is x^2, and the integer case is the dangerous one: 3^2 returns 1 silently, with no error to warn you.

Mistake 2: Mixing The Math Module With SymPy Symbols

Where it slips in: Importing math out of habit, then building a SymPy expression with math.sin(x).

Don't do this: Mixing libraries. math.sin(x) with a symbolic x stops with TypeError: Cannot convert expression to float, because math only accepts numbers.

The correct way: Use sympy.sin, sympy.exp and sympy.log for expressions, and np.sin only on arrays after lambdify.

Mistake 3: Copying SymPy's Antiderivative Without + C

Where it slips in: Writing up an indefinite integral straight from integrate(x**2, x).

Don't do this: Reporting $\int x^2,dx = \frac{x^3}{3}$. SymPy returns one antiderivative and omits the constant on purpose.

The correct way: Write $\frac{x^3}{3} + C$, the family of all antiderivatives, as covered in constant of integration.

The memorizer, who learned "integrate means add C", gets this right on paper and then drops the constant the moment a computer prints the answer. The output looks finished, so it gets trusted.

Mistake 4: Shrinking The Step Size Until The Numbers Break

Where it slips in: A home-made numerical derivative where a smaller step "must" be more accurate.

Don't do this: Setting h = 1e-14. The step-size table above shows the error jumping back up to $3.71 \times 10^{-3}$.

The correct way: Use the central difference with $h$ around $10^{-5}$, or skip approximation entirely with diff plus lambdify. Check for a zero step before dividing, not after; a reply on that same Physics Forums thread caught code that tested whether h was 0 only after it had already been used.

The rusher reaches for a tiny h because the limit definition says $h \to 0$. Students who first print an error table, like the one above, stop trusting "smaller is better" and start choosing $h$ on purpose.

Practice Problems On Calculus With Python

Solve each by hand first, then confirm in Python.

1. Use limit to find $\lim_{x \to 3} \frac{x^2 - 9}{x - 3}$.

2. Use diff to differentiate $\ln\left(x^2 + 1\right)$.

3. Use integrate to evaluate $\int_{0}^{\pi} \sin x,dx$.

4. Use series to find the Maclaurin polynomial of $\cos x$ up to the $x^4$ term.

5. Use dsolve to solve $\frac{dy}{dt} = 3y$ with $y(0) = 2$.

6. Use np.trapezoid with np.linspace(0, 1, 11) to approximate $\int_{0}^{1} x^2,dx$, and compare with the exact value.

Answers (each run in SymPy 1.14 or NumPy 2.5):

• Answer to Question 1: $6$.

• Answer to Question 2: $\frac{2x}{x^2 + 1}$.

• Answer to Question 3: $2$.

• Answer to Question 4: $1 - \frac{x^2}{2} + \frac{x^4}{24}$.

• Answer to Question 5: $y = 2e^{3t}$.

• Answer to Question 6: $0.335$ (printed as 0.33499999999999996), against the exact $\frac{1}{3} \approx 0.3333$.

Wrapping Up

• Calculus with Python splits into symbolic work in SymPy, which returns exact formulas, and numerical work in NumPy or SciPy, which returns fast decimals.

• limit, diff, integrate, series and dsolve cover limits, derivatives, integrals, Taylor series and differential equations in one call each.

• SymPy omits $+ C$ on indefinite integrals, so the constant is yours to add.

• Floats fail quietly: a tiny step $h$ or a huge $n$ can turn a correct method into a wrong number.

• Write powers as **, never ^, and use SymPy's own functions on symbols.

Your Next Move

Work through problems 1 to 3 by hand first, then confirm each with one line of SymPy; if problem 6 gives a number you did not expect, return to the section on integrating numerically with NumPy. To practise these with a live Bhanzu trainer, book a free demo class

Read More

• Calculus formulas reference

• Table of derivatives

• Table of integrals

• Maclaurin series

• Second derivative test