Taupy, Adventures in Math with Symbolic Programming

Formulating a method as a computer-executable program and debugging that program is a powerful exercise in the learning process.

— gjs

Introduction

Math definitions concisely expressed with SymPy, a symbolic Python library.

Calculus

Sequence

Adjacent Sequence

Two adjacent sequences converge and have the same limit.

from sympy import *

n = symbols("n")

un = 1 - 1 / n
vn = 1 + 1 / n

is_adjacent = limit(un - vn, n, oo).doit() == 0
assert is_adjacent

is_same_limit = limit(un, n, oo).doit() == limit(vn, n, oo).doit()
assert is_same_limit

Cauchy Sequence

A set is complete if and only if every Cauchy sequence converges within it, hence Q is not complete.

from sympy import *

n = symbols("n")

cauchy_sequence_limit = summation(1 / factorial(n), (n, 0, oo))

is_limit_not_in_Q = cauchy_sequence_limit == E
assert is_limit_not_in_Q

Cesaro Mean

When a sequence has a limit the Cesaro mean has the same limit.

from sympy import *

n = symbols("n")

expression = n * log(1 + 2 / n)
sequence_limit = limit(expression, n, oo).doit()

cesaro_mean = 0
N = 1000
for i in range(1, N):
    cesaro_mean += expression.subs({n: i})
cesaro_mean /= N

is_same_limit = abs(cesaro_mean.evalf() - sequence_limit) < 0.02
assert is_same_limit

Recursive Sequence

The set of solutions to a second-order constant-recursive sequence is a vector space.

from sympy import *

r, un1, un, n, c, d = symbols("r un1 un n c d")

a = -4
b = -3

un2 = a * un1 + b * un
characteristic_equation_of_un2 = r**2 - a * r - b
solutions = solve(characteristic_equation_of_un2, r)

set_of_solutions = c * solutions[0] ** n + d * solutions[1] ** n

u0 = set_of_solutions.subs({c: 50, d: 100, n: 0})
u1 = set_of_solutions.subs({c: 50, d: 100, n: 1})
u2 = un2.subs({un1: u1, un: u0})
is_solution = u2 == set_of_solutions.subs({c: 50, d: 100, n: 2})
assert is_solution

mu = 2
u10 = set_of_solutions.subs({c: 50, d: 100, n: 10})
u10_prime = set_of_solutions.subs({c: -50, d: -100, n: 10})
v10 = set_of_solutions.subs({c: 60, d: 110, n: 10})
w10 = set_of_solutions.subs({c: 70, d: 150, n: 10})
neutral_zero = set_of_solutions.subs({c: 0, d: 0, n: 10})
neutral_one = set_of_solutions.subs({c: 0, d: 1, n: 10})

is_associative = (
    (u10 + v10) + w10
    == u10 + (v10 + w10)
    == set_of_solutions.subs({c: 50 + 60 + 70, d: 100 + 110 + 150, n: 10})
)
is_commutative = (
    (u10 + v10)
    == (v10 + u10)
    == set_of_solutions.subs({c: 50 + 60, d: 100 + 110, n: 10})
)
is_distributive = (
    mu * (u10 + v10)
    == mu * u10 + mu * v10
    == set_of_solutions.subs({c: mu * (50 + 60), d: mu * (100 + 110), n: 10})
)
is_symmetry = u10 + u10_prime == neutral_zero
is_neutral_add = u10 + neutral_zero == u10
is_neutral_mult = u10 * neutral_one == u10

is_vector_space = (
    is_associative
    and is_commutative
    and is_distributive
    and is_symmetry
    and is_neutral_add
    and is_neutral_mult
)
assert is_vector_space

Integral

Signal Value

The mean value of a signal is the integral of its amplitude over time divided by the duration.

from sympy import *

t = symbols("t")

f = sin(t) ** 2

t1 = 0
t2 = 2 * pi

mean_value = integrate(f, (t, t1, t2)) / (t2 - t1)

is_compliant_with_plot = abs(mean_value - 0.5) < 0.001
assert is_compliant_with_plot

Convolution Operator

The convolution operator combines two functions by integrating the product of one function with a shifted and reversed version of the other.

from sympy import *

x, t = symbols("x t")

rectangular_function_f = Piecewise((1, And(t >= -1 / 2, t <= 1 / 2)), (0, True))
shifted_rectangular_function_g = Piecewise(
    (1, And(x - t >= -1 / 2, x - t <= 1 / 2)), (0, True)
)
convolution_f_g = integrate(
    rectangular_function_f * shifted_rectangular_function_g, (t, -oo, oo)
)

is_compliant_with_plot = all(
    [
        convolution_f_g.subs({x: -0.25}) == 0.75,
        convolution_f_g.subs({x: 0}) == 1.0,
        convolution_f_g.subs({x: 0.25}) == 0.75,
        convolution_f_g.subs({x: 0.5}) == 0.5,
        convolution_f_g.subs({x: 1}) == 0,
    ]
)
assert is_compliant_with_plot

Cauchy Principal

The Cauchy principal value assigns values to certain improper integrals that would otherwise be undefined.

from sympy import *

x = symbols("x")
r, epsilon = symbols("r epsilon", positive=True)

f = 1 / x

is_undefined = integrate(f, (x, -oo, oo)) == nan
assert is_undefined

undefined_at_x = 0

cauchy_principal_value = simplify(
    integrate(f, (x, -r, undefined_at_x - epsilon))
    + integrate(f, (x, undefined_at_x + epsilon, r))
)

is_defined = cauchy_principal_value != nan
assert is_defined

Jacobian Determinant

The determinant of the Jacobian matrix provides the scaling factor that adjusts the integral to the new coordinate system.

from sympy import *

x, y, r, theta = symbols("x y r theta")

xp = r * cos(theta)
yp = r * sin(theta)

jacobian_matrix = Matrix(
    [[diff(xp, r), diff(xp, theta)], [diff(yp, r), diff(yp, theta)]]
)
jacobian_polar_factor = det(jacobian_matrix).simplify()

integrate_fxy_in_unit_circle_region_polar = integrate(
    (xp - yp) ** 2 * jacobian_polar_factor, (r, 0, 1), (theta, 0, 2 * pi)
)

integrate_fxy_in_unit_circle_region_cartesian = integrate(
    (x - y) ** 2, (y, -sqrt(1 - x**2), sqrt(1 - x**2)), (x, -1, 1)
)

is_same_area = (
    integrate_fxy_in_unit_circle_region_polar
    == integrate_fxy_in_unit_circle_region_cartesian
    == pi / 2
)
assert is_same_area

Approximation & Interpolation

Newton's Method

Newton's method is an iterative technique for finding the roots of a function by using successive approximations.

from sympy import *

x = symbols("x")

f = x**2 - 4
Df = diff(f, x)

x0 = 10
x1 = x0 - f.subs({x: x0}) / Df.subs({x: x0})
x2 = x1 - f.subs({x: x1}) / Df.subs({x: x1})
x3 = x2 - f.subs({x: x2}) / Df.subs({x: x2})

tangent_at_x0 = f.subs({x: x0}) + Df.subs({x: x0}) * (x - x0)

is_compliant_with_plot = all(
    [
        x0 == 10,
        abs(x1 - 5.2) < 0.01,
        abs(x2 - 2.984) < 0.01,
        abs(x3 - 2.162) < 0.01,
        tangent_at_x0.subs({x: x1}) == 0,
    ]
)
assert is_compliant_with_plot

Ordinary Least Squares

Ordinary Least Squares finds the best-fitting line by minimizing prediction errors across data points.

from sympy import *

x1, x2 = symbols("x1 x2")
c1 = 3
c2 = 2
c3 = 5

f = c1 + c2 * x1 + c3 * x2

X = Matrix([[1, 0, 1], [1, 1, 2], [1, 2, 3], [1, 3, 5], [1, 4, 7]])
y = Matrix(
    [
        f.subs({x1: 0, x2: 1}),
        f.subs({x1: 1, x2: 2}),
        f.subs({x1: 2, x2: 3}),
        f.subs({x1: 3, x2: 5}),
        f.subs({x1: 4, x2: 7}),
    ]
)

beta_hat = (X.T * X).inv() * X.T * y
y_hat = X * beta_hat

is_good_estimator = beta_hat == Matrix([[c1], [c2], [c3]]) and y_hat == y
assert is_good_estimator

Lagrange Polynomial

Lagrange polynomial interpolates a set of points by combining simpler polynomial terms.

from sympy import *

x = symbols("x")

f = x**3 - 2 * x - 5

x0 = 0.5
y0 = f.subs({x: x0})
x1 = 2.0
y1 = f.subs({x: x1})
x2 = 3.5
y2 = f.subs({x: x2})

lagrange_basis = [
    ((x - x1) * (x - x2)) / ((x0 - x1) * (x0 - x2)),
    ((x - x0) * (x - x2)) / ((x1 - x0) * (x1 - x2)),
    ((x - x0) * (x - x1)) / ((x2 - x0) * (x2 - x1)),
]

lagrange_polynomial = (
    y0 * lagrange_basis[0] + y1 * lagrange_basis[1] + y2 * lagrange_basis[2]
)

is_compliant_with_plot = all(
    [
        abs(lagrange_polynomial.subs({x: x0}) + 5.875) < 0.01,
        abs(lagrange_polynomial.subs({x: x1}) + 1.0) < 0.01,
        abs(lagrange_polynomial.subs({x: x2}) - 30.875) < 0.01,
    ]
)
assert is_compliant_with_plot

Chebyshev Nodes

Chebyshev nodes are the roots of Chebyshev polynomials used for optimal polynomial interpolation.

from sympy import *

x, k, n, a, b = symbols("x k n a b")

f = x**3 - 2 * x - 5

theta = acos(x)
chebyshev_polynomial = cos(n * theta)

chebyshev_polynomial_definition = (
    cos(0) == 1
    and cos(theta) == x
    and simplify(cos(2 * theta) - (2 * cos(theta) ** 2 - 1)) == 0
    and 2 * cos(theta) ** 2 - 1
    == 2 * x**2 - 1
    == 2 * x * chebyshev_polynomial.subs({n: 1}) - chebyshev_polynomial.subs({n: 0})
    and integrate(
        chebyshev_polynomial.subs({n: 1})
        * chebyshev_polynomial.subs({n: 2})
        / sqrt(1 - x**2),
        (x, -1, 1),
    )
    == 0
)
assert chebyshev_polynomial_definition

chebyshev_node = (a + b) / 2 + ((b - a) / 2) * cos((2 * k + 1) / (2 * n) * pi)
chebyshev_node3 = [
    chebyshev_node.subs({k: 0, n: 3, a: -1, b: 1}),
    chebyshev_node.subs({k: 1, n: 3, a: -1, b: 1}),
    chebyshev_node.subs({k: 2, n: 3, a: -1, b: 1}),
]

t3 = chebyshev_polynomial.subs({n: 3})
is_t3_roots = (
    t3.subs({x: chebyshev_node3[0]})
    == t3.subs({x: chebyshev_node3[1]})
    == t3.subs({x: chebyshev_node3[2]})
    == 0
)
assert is_t3_roots

x0 = chebyshev_node.subs({k: 0, n: 3, a: 0, b: 4})
y0 = f.subs({x: x0})
x1 = chebyshev_node.subs({k: 1, n: 3, a: 0, b: 4})
y1 = f.subs({x: x1})
x2