Python (SymPy) Programming for solving equations and number sequences

Here I'd like to share how to solve equations using Python, in particular "SymPy", a Python library for symbolic formula manipulation.

In addition to (simultaneous) equations, I'd like to show you how to find a number sequence defined by a recurrence relation. The following article provides an essence and a flavor of how to use SymPy:

pianofisica.hatenablog.com

Algebraic equation
Numerical solutions
Simultaneous linear equations
- 2 variables
- 3 variables
Simultaneous equations
Recurrence relation

Algebraic equation

An algebraic equation is a relation equated between polynomials, which is made up of finite numbers of additions and multiplications of variables and constants.

A number which solves an algebraic equation is called "algebraic number". Thus, integers, rational numbers, irrational numbers like $\sqrt{2}$ are all algebraic numbers.

Whereas those are not algebraic numbers are called "transcendental numbers". Correspondingly, equations that are not algebraic are called "transcendental equations".

Examples for transcendental equations are

$1-\log x=0 \qquad \sin x=0$

Solutions to these equations, $e$ Napier's constant and $\pi$ pi, are thus "transcendental numbers". In the following, I'd like to show you how to solve algebraic equations using SymPy.

Linear equation (1st-degree polynomial equation)

To solve the equation

$ax+b=0$

input the following Python source

import sympy as sp
sp.var('x, a, b')
Sol1=sp.solve (a*x+b, x)
sp.init_printing()
display(Sol1[0])

then you get the answer as its output

$\displaystyle{x=-\frac{b}{a}}$

Quadratic equation (2nd-degree polynomial equation)

To solve the equation

$ax^2+bx+c=0$

input the following Python source

sp.var('x, a, b, c')
Sol2=sp.solve (a*x**2+b*x+c, x)
display(Sol2)

then you get the answer as its output

$\displaystyle{x={-b-{\sqrt{b^2-4ac}}\over{2a}} ,\qquad x={-b+{\sqrt{b^2-4ac}}\over{2a}}}$

Cubic equation (3rd-degree polynomial equation)

To solve the equation

$ax^3+bx^2+cx+d=0$

input the following Python source

sp.var('x, a, b, c, d')
Sol3=sp.solve (a*x**3+b*x**2+c*x+d, x)

display(Sol3[0])

display(Sol3[1])

display(Sol3[2])

then you get the answer as its output

$\displaystyle{x=r^{1/3}-\frac{s}{r^{1/3}}-\frac{b}{3a},\qquad x=\left(\frac{-1\mp i\sqrt{3}}{2}\right)r^{1/3}-\left(\frac{-1\pm i\sqrt{3}}{2}\right)\frac{s}{r^{1/3}}-\frac{b}{3a}\\ r={{\left(27a^2d^2+\left(4b^3-18abc\right)d+4ac^3-b^2c^2\right)^{1/2}}\over{2\,3^{3/2}a^2}}+{{{{bc}}-{{3ad}}}\over{6a^2}}-{{b^3}\over{27a^3}} \\ s=-\frac{b^2}{9a^2}+\frac{c}{3a} }$

This is the formula for roots of the cubic equation.

Quartic equation (4th-degree polynomial equation)

To solve the equation

$ax^4+bx^3+cx^2+dx+e=0$

input the following Python source

sp.var('x, a, b, c, d, e')
Sol4=sp.solve (a*x**4+b*x**3+c*x**2+d*x+e, x)
display(Sol4)

then you get the answer as its output,
but the answer has lengthy expression, I'd omit here the result.

Quintic equation (5th-degree polynomial equation)

$ax^4+bx^3+cx^2+dx+e=0$

is known to be never solved by algebraic manipulations.
Indeed, even if you input

sp.var('x, a, b, c, d, e, f')
Sol5=sp.solve (a*x**5+b*x**4+c*x**3+d*x**2+e*x+f, x)
display(Sol5)

there's not output.

Numerical solutions

So far we have seen how to solve an algebraic equation for a variable $\ x$ ,
in general, no equation of order more than 5 can be solved algebraically.

However, for some purpose, it is sometimes enough to know a root numerically:

For example, the equation

$x^5+6x^4+x^3+x+1=0$

is of order 5 and thus cannot be solved in algebraic way:

sp.var('x')
sp.solve (x**5+6*x**4+x**3+x+1, x)

Whereas, using "NumPy" a Python library for numerical calculation,

import numpy
numpy.roots([1, 6, 1, 0, 1, 1])

one can find its roots numerically

$x_1=-5.8241071 \\ x_2 =0.42084104+0.52540961\, i \\ x_3= 0.42084104-0.52540961\,i \\ x_4=-0.50878749+0.34645119\,i \\x_5=-0.50878749-0.34645119\,i$

Simultaneous linear equations

2 variables

Simultaneous equations for two variables $x$ and $y$

$ax+by=e \\ cx+dy=f$

can be solved by Python via

sp.var('x, y, a, b, c, d, e, f')
eq1=sp.Eq(a*x+b*y, e)
eq2=sp.Eq(c*x+d*y, f)
sp.solve ([eq1, eq2], [x, y])

then you could get the answer

$\displaystyle{x={{de-bf}\over{ad-bc}} ,\qquad y={{af-ce}\over{ad-bc}}}$

for example

$2x+y=3 \\ x+3y=4$

is solved by

sp.var('x, y')
eq3=sp.Eq(2*x+1*y, 3)
eq4=sp.Eq(1*x+3*y, 4)
sp.solve ([eq3, eq4], [x, y])

and thus

$x=1, \quad y=1$

3 variables

Similarly, for the case with three variables $x$ , $y$ and $z$ :

$ax+by+cz=j \\ dx+ey+fz=k \\ gx+hy+iz=l$

is solved by

sp.var('x, y, z, a, b, c, d, e, f, g, h, i, j, k, l')
eq1=sp.Eq(a*x+b*y+c*z, j)
eq2=sp.Eq(d*x+e*y+f*z, k)
eq3=sp.Eq(g*x+h*y+i*z, l)
sp.solve ([eq1, eq2, eq3], [x, y, z])

and thus

$\displaystyle{x={{j(ei-fh)-k(bi-ch)+l(bf-ce)}\over{a(ei-fh)-b(di-fg)+c(dh-eg)}} \\ y={{−𝑗(𝑑𝑖−𝑓𝑔)+𝑘(𝑎𝑖−𝑐𝑔)−𝑙(𝑎𝑓−𝑐𝑑)}\over{a(ei-fh)-b(di-fg)+c(dh-eg)}} \\ z={{𝑗(𝑑ℎ−𝑒𝑔)−𝑘(𝑎ℎ−𝑏𝑔)+𝑙(𝑎𝑒−𝑏𝑑)}\over{a(ei-fh)-b(di-fg)+c(dh-eg)} }}$

Simultaneous linear equations can also be represented in terms of matrix,
and using its inverse matrix, one can also solve the equations.

Simultaneous equations

For general simultaneous equation can also be solved by

# "sympy.solve([eq1, eq2, ...], [var1, var2, ...])" solves eq1, eq2, ...

For example,
a cross point of parabola $y=ax^2+bx+c$ and a line $y=gx+h$ is given by a solution to

$y=ax^2+bx+c \\ y=gx+h$

which is solved by Python via

sp.var('x, y, a, b, c, g, h')
eq5=sp.Eq(y, a*x**2+b*x+c)
eq6=sp.Eq(y, g*x+h)
sp.solve ([eq5, eq6], [x, y])

Recurrence relation

A number sequence defined by a recurrence relation like

$F_{n+2}=F_{n+1}+F_n \qquad F_0=0 \qquad F_1 =1$

can be found by using Python as follows:

def n(k):
    if k == 0:
        return 0
    if k == 1:
        return 1
    else:
        return n(k-1) + n(k-2) 
N_list = []
for i in range(10):
    N = n(i)
    N_list.append(N)
print(N_list)