Initialize Enviroment

In [1]:

from IPython.display import display, Math, Latex 
from sympy import *
init_printing()

from helper import comparator_factory, comparator_eval_factory, comparator_method_factory

x,y,z = symbols('x y z')

func_comparator = comparator_factory('Before applying {}:','After:')
method_comparator = comparator_method_factory('Before calling {}:','After:')
eval_comparator = comparator_eval_factory('Before evaluation:','After evaluation:')

Matrices

Construct matrix

Construct matrix from list

To make a matrix in Sympy, initialize a Matrix class by providing a list of row vectors (nested list).

In [2]:

m = Matrix(
    [
        [1, -1], 
        [3, 4], 
        [0, 2]
    ]
)

m

Out[2]:

$\displaystyle \left[\begin{matrix}1 & -1\\3 & 4\\0 & 2\end{matrix}\right]$

Pass a single layer list of elements will create a $n \times 1$ matrix.

In [3]:

m = Matrix([1,2,3])

m

Out[3]:

$\displaystyle \left[\begin{matrix}1\\2\\3\end{matrix}\right]$

Common Matrix Constructors

Several constructors exist for creating common matrices.

Identity Matrix

To create identity matrix, use eye(n)

In [4]:

M = eye(3)

M

Out[4]:

$\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]$

Matrix of Zero

To create a zero matrix, use zero(n)

In [5]:

M = zeros(3,4)

display(M)

$\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\0 & 0 & 0 & 0\\0 & 0 & 0 & 0\end{matrix}\right]$

Matrix of Ones

To create a zero matrix, use one(n)

In [6]:

M = ones(3,2)

M

Out[6]:

$\displaystyle \left[\begin{matrix}1 & 1\\1 & 1\\1 & 1\end{matrix}\right]$

Diagonal

Pass multiple matrix or numbers to construct a diagonal matrix. Numbers will be intepreted as $1 \times 1$ matrix.

In [7]:

M = diag(-1, ones(2, 2), Matrix([5, 7, 5]))

M

Out[7]:

$\displaystyle \left[\begin{matrix}-1 & 0 & 0 & 0\\0 & 1 & 1 & 0\\0 & 1 & 1 & 0\\0 & 0 & 0 & 5\\0 & 0 & 0 & 7\\0 & 0 & 0 & 5\end{matrix}\right]$

Basic operations

Shape

To get the shape of a matrix, use shape

In [8]:

M.shape

Out[8]:

$\displaystyle \left( 6, \ 4\right)$

Get rows and columns

Use row() and col() to get single row and column.

To get the first row

In [9]:

M.row(0)

Out[9]:

$\displaystyle \left[\begin{matrix}-1 & 0 & 0 & 0\end{matrix}\right]$

To get the last column

In [10]:

M.col(-1)

Out[10]:

$\displaystyle \left[\begin{matrix}0\\0\\0\\5\\7\\5\end{matrix}\right]$

Delete Rows and Columns

To delete a row or column, use row_del and col_del. These operaitons modify the matrix inplace.

For example, we create a matrix then delete the first row and last column.

In [11]:

M = Matrix([[1, 2, 3], [3, 2, 1]])

display('Before deletion:',M)

M.row_del(0)
M.col_del(-1)

display('After deletion:',M)

'Before deletion:'

$\displaystyle \left[\begin{matrix}1 & 2 & 3\\3 & 2 & 1\end{matrix}\right]$

'After deletion:'

$\displaystyle \left[\begin{matrix}3 & 2\end{matrix}\right]$

Transpose of Matrix

Use T

In [12]:

M = Matrix([[1, 2, 3], [4, 5, 6]])

display('Original matrix:',M)
display('Transposed matrix:',M.T)

'Original matrix:'

$\displaystyle \left[\begin{matrix}1 & 2 & 3\\4 & 5 & 6\end{matrix}\right]$

'Transposed matrix:'

$\displaystyle \left[\begin{matrix}1 & 4\\2 & 5\\3 & 6\end{matrix}\right]$

Insert Rows and Columns

To insert rows and columns, use row_insert and col_insert. These operations don't operate inplace and return a new matrix with modified data.

row_insert() accepts a position and a $1 \times n$ matrix.

In [13]:

M = Matrix([[1, 2, 3], [3, 2, 1]])

display('Before insertion:',M)

M = M.row_insert(0, Matrix([[6,6,6]]))

display('After Insertion:',M)

'Before insertion:'

$\displaystyle \left[\begin{matrix}1 & 2 & 3\\3 & 2 & 1\end{matrix}\right]$

'After Insertion:'

$\displaystyle \left[\begin{matrix}6 & 6 & 6\\1 & 2 & 3\\3 & 2 & 1\end{matrix}\right]$

col_insert() accepts position and $n \times 1$ matrix.

In [14]:

M = Matrix([[1, 2, 3], [3, 2, 1]])

display('Before insertion:',M)

M = M.col_insert(-1, Matrix([6,6]))

display('After Insertion:',M)

'Before insertion:'

$\displaystyle \left[\begin{matrix}1 & 2 & 3\\3 & 2 & 1\end{matrix}\right]$

'After Insertion:'

$\displaystyle \left[\begin{matrix}1 & 2 & 6 & 3\\3 & 2 & 6 & 1\end{matrix}\right]$

Matrix Computation

Basic matrix computation like addition, substraction, multiplication and inversion can be achived through python operators +,-,*,**

Addition

In [15]:

M = Matrix([[1, 3], [-2, 3]])
N = Matrix([[0, 3], [0, 7]])

print('M')
display(M)
print('N')
display(N)
print('M+N')
display(M+N)

$\displaystyle \left[\begin{matrix}1 & 3\\-2 & 3\end{matrix}\right]$

$\displaystyle \left[\begin{matrix}0 & 3\\0 & 7\end{matrix}\right]$

M+N

$\displaystyle \left[\begin{matrix}1 & 6\\-2 & 10\end{matrix}\right]$

Substraction

In [16]:

M = Matrix([[1, 3], [-2, 3]])
N = Matrix([[0, 3], [0, 7]])

print('M')
display(M)
print('N')
display(N)
print('M-N')
display(M-N)

$\displaystyle \left[\begin{matrix}1 & 3\\-2 & 3\end{matrix}\right]$

$\displaystyle \left[\begin{matrix}0 & 3\\0 & 7\end{matrix}\right]$

M-N

$\displaystyle \left[\begin{matrix}1 & 0\\-2 & -4\end{matrix}\right]$

Multiplication

In [17]:

M = Matrix([[1, 2, 3], [3, 2, 1]])
N = Matrix([0, 1, 1])

print('M')
display(M)
print('N')
display(N)
display('MxN')
display(M*N)

$\displaystyle \left[\begin{matrix}1 & 2 & 3\\3 & 2 & 1\end{matrix}\right]$

$\displaystyle \left[\begin{matrix}0\\1\\1\end{matrix}\right]$

'MxN'

$\displaystyle \left[\begin{matrix}5\\3\end{matrix}\right]$

Inverse

In [18]:

M = Matrix([[1, 3], [-2, 3]])

print('M')
display(M)
print('inverse of M:')
display(M**(-1))

$\displaystyle \left[\begin{matrix}1 & 3\\-2 & 3\end{matrix}\right]$

inverse of M:

$\displaystyle \left[\begin{matrix}\frac{1}{3} & - \frac{1}{3}\\\frac{2}{9} & \frac{1}{9}\end{matrix}\right]$

Determinant

To get determinant of matrix, use det

In [19]:

M = Matrix([[1, 0, 1], [2, -1, 3], [4, 3, 2]])

method_comparator(M,'det')

Before calling det():

$\displaystyle \left[\begin{matrix}1 & 0 & 1\\2 & -1 & 3\\4 & 3 & 2\end{matrix}\right]$

After:

$\displaystyle -1$

RREF

To put a matrix into reduced echelon form, use rref(). It returns a matrix of two elements, the first is the reduced echelon form and the second is a list of indices of the pivot columns.

In [20]:

M = Matrix([[1, 0, 1, 3], [2, 3, 4, 7], [-1, -3, -3, -4]])

M_rref, pivot_columns = M.rref()

print('Original Matrix:')
display(M)

print('Reduced echelon form:')
display(M_rref)

print('Pivot columns:')
display(pivot_columns)

Original Matrix:

$\displaystyle \left[\begin{matrix}1 & 0 & 1 & 3\\2 & 3 & 4 & 7\\-1 & -3 & -3 & -4\end{matrix}\right]$

Reduced echelon form:

$\displaystyle \left[\begin{matrix}1 & 0 & 1 & 3\\0 & 1 & \frac{2}{3} & \frac{1}{3}\\0 & 0 & 0 & 0\end{matrix}\right]$

Pivot columns:

$\displaystyle \left( 0, \ 1\right)$

Nullspace

nullspace() returns a list of column vectors which span the nullspace of a matrix.

In [21]:

M = Matrix([[1, 2, 3, 0, 0], [4, 10, 0, 0, 1]])

print('Matrix M:')
display(M)
print('Spanning vectors for nullspace of M:')
display(M.nullspace())

Matrix M:

$\displaystyle \left[\begin{matrix}1 & 2 & 3 & 0 & 0\\4 & 10 & 0 & 0 & 1\end{matrix}\right]$

Spanning vectors for nullspace of M:

$\displaystyle \left[ \left[\begin{matrix}-15\\6\\1\\0\\0\end{matrix}\right], \ \left[\begin{matrix}0\\0\\0\\1\\0\end{matrix}\right], \ \left[\begin{matrix}1\\- \frac{1}{2}\\0\\0\\1\end{matrix}\right]\right]$

Column space

columnspace() returns a list of column vectors which span the columnspace of a matrix.

In [22]:

M = Matrix([[1, 1, 2], [2 ,1 , 3], [3 , 1, 4]])

print('Matrix M:')
display(M)
print('Spanning vectors for columnspace of M:')
display(M.columnspace())

Matrix M:

$\displaystyle \left[\begin{matrix}1 & 1 & 2\\2 & 1 & 3\\3 & 1 & 4\end{matrix}\right]$

Spanning vectors for columnspace of M:

$\displaystyle \left[ \left[\begin{matrix}1\\2\\3\end{matrix}\right], \ \left[\begin{matrix}1\\1\\1\end{matrix}\right]\right]$

Eigenvalues and Eigenvectors

To find the eigenvalues of a matrix, use eigenvals. It returns a dictionary of eigenvalue:algebraic multiplicity pairs.

In [23]:

M = Matrix([[3, -2,  4, -2], [5,  3, -3, -2], [5, -2,  2, -2], [5, -2, -3,  3]])

print('Matrix M')
display(M)

print('Eigenvalues and their algebraic multiplicity:')
display(M.eigenvals())

Matrix M

$\displaystyle \left[\begin{matrix}3 & -2 & 4 & -2\\5 & 3 & -3 & -2\\5 & -2 & 2 & -2\\5 & -2 & -3 & 3\end{matrix}\right]$

Eigenvalues and their algebraic multiplicity:

$\displaystyle \left\{ -2 : 1, \ 3 : 1, \ 5 : 2\right\}$

eigenvectors return a list of tuples of the form (eigenvalue, algebraic multiplicity, [eigenvectors])

In [24]:

M = Matrix(
    [
        [3, -2,  4, -2], 
        [5,  3, -3, -2], 
        [5, -2,  2, -2], 
        [5, -2, -3,  3]
    ]
)

print('Matrix M')
display(M)

print('Eigenvalues, their algebraic multiplicity and eigenvectors:')
display(M.eigenvects())

Matrix M

$\displaystyle \left[\begin{matrix}3 & -2 & 4 & -2\\5 & 3 & -3 & -2\\5 & -2 & 2 & -2\\5 & -2 & -3 & 3\end{matrix}\right]$

Eigenvalues, their algebraic multiplicity and eigenvectors:

$\displaystyle \left[ \left( -2, \ 1, \ \left[ \left[\begin{matrix}0\\1\\1\\1\end{matrix}\right]\right]\right), \ \left( 3, \ 1, \ \left[ \left[\begin{matrix}1\\1\\1\\1\end{matrix}\right]\right]\right), \ \left( 5, \ 2, \ \left[ \left[\begin{matrix}1\\1\\1\\0\end{matrix}\right], \ \left[\begin{matrix}0\\-1\\0\\1\end{matrix}\right]\right]\right)\right]$

If all you want is the characteristic polymonial, use charpoly() and factor()

In [25]:

print('Matrix M')
display(M)

lamda = symbols('lamda')
p = M.charpoly(lamda).factor()

print('characteristic polymonial of M')
display(p)

Matrix M

$\displaystyle \left[\begin{matrix}3 & -2 & 4 & -2\\5 & 3 & -3 & -2\\5 & -2 & 2 & -2\\5 & -2 & -3 & 3\end{matrix}\right]$

characteristic polymonial of M

$\displaystyle \left(\lambda - 5\right)^{2} \left(\lambda - 3\right) \left(\lambda + 2\right)$

Diagonalization

To diagonalize a matrix , use diagonalize(). diagonalize returns a tuple ($P$,$D$), where D is diagonal and $M = PDP^{-1}$.

In [26]:

M = Matrix(
    [
        [3, -2,  4, -2], 
        [5,  3, -3, -2], 
        [5, -2,  2, -2], 
        [5, -2, -3,  3]
    ]
)

P, D = M.diagonalize()

print('Matrix M')
display(M)

print('Matrix P')
display(P)

print('Matrix D')
display(D)

display(Latex('$PDP^{-1}$'))
display(P*D*P**-1)

Matrix M

$\displaystyle \left[\begin{matrix}3 & -2 & 4 & -2\\5 & 3 & -3 & -2\\5 & -2 & 2 & -2\\5 & -2 & -3 & 3\end{matrix}\right]$

Matrix P

$\displaystyle \left[\begin{matrix}0 & 1 & 1 & 0\\1 & 1 & 1 & -1\\1 & 1 & 1 & 0\\1 & 1 & 0 & 1\end{matrix}\right]$

Matrix D

$\displaystyle \left[\begin{matrix}-2 & 0 & 0 & 0\\0 & 3 & 0 & 0\\0 & 0 & 5 & 0\\0 & 0 & 0 & 5\end{matrix}\right]$

$PDP^{-1}$

$\displaystyle \left[\begin{matrix}3 & -2 & 4 & -2\\5 & 3 & -3 & -2\\5 & -2 & 2 & -2\\5 & -2 & -3 & 3\end{matrix}\right]$

Solvers

Basic equations

The main function for solving algebraic equations is solveset(). The syntax for solveset() is solveset(equation, variable=None, domain=S.Complexes).

In [27]:

eq = Eq(x**2-1, 0)

func_comparator(eq, solveset, x)

Before applying solveset():

$\displaystyle x^{2} - 1 = 0$

After:

$\displaystyle \left\{-1, 1\right\}$

If an expresison instead of a Eq instance is passed in solveset(), it is automatically assumed to be zero.

In [28]:

expr = x**2-1

func_comparator(eq, solveset, x)

Before applying solveset():

$\displaystyle x^{2} - 1 = 0$

After:

$\displaystyle \left\{-1, 1\right\}$

Result of the solver can be a infinite set.

In [29]:

eq = Eq(sin(x) - 1)

func_comparator(eq, solveset, x)

Before applying solveset():

/usr/local/lib/python3.7/site-packages/sympy/core/relational.py:470: SymPyDeprecationWarning: 

Eq(expr) with rhs default to 0 has been deprecated since SymPy 1.5.
Use Eq(expr, 0) instead. See
https://github.com/sympy/sympy/issues/16587 for more info.

  deprecated_since_version="1.5"

$\displaystyle \sin{\left(x \right)} - 1 = 0$

After:

$\displaystyle \left\{2 n \pi + \frac{\pi}{2}\; |\; n \in \mathbb{Z}\right\}$

Result of the solver can be an empty set.

In [30]:

eq = Eq(exp(x))

func_comparator(eq, solveset, x)

Before applying solveset():

/usr/local/lib/python3.7/site-packages/sympy/core/relational.py:470: SymPyDeprecationWarning: 

Eq(expr) with rhs default to 0 has been deprecated since SymPy 1.5.
Use Eq(expr, 0) instead. See
https://github.com/sympy/sympy/issues/16587 for more info.

  deprecated_since_version="1.5"

$\displaystyle \text{False}$

After:

$\displaystyle \emptyset$

When solver fails to find solutions, it returns an conditional set.

In [31]:

eq = Eq(cos(x) - x)

func_comparator(eq, solveset, x)

Before applying solveset():

/usr/local/lib/python3.7/site-packages/sympy/core/relational.py:470: SymPyDeprecationWarning: 

Eq(expr) with rhs default to 0 has been deprecated since SymPy 1.5.
Use Eq(expr, 0) instead. See
https://github.com/sympy/sympy/issues/16587 for more info.

  deprecated_since_version="1.5"

$\displaystyle - x + \cos{\left(x \right)} = 0$

After:

$\displaystyle \left\{x \mid x \in \mathbb{C} \wedge - x + \cos{\left(x \right)} = 0 \right\}$

solveset() reports each solution only once.

In [32]:

eq = Eq(x**3 - 6*x**2 + 9*x,0)

func_comparator(eq, solveset, x)

Before applying solveset():

$\displaystyle x^{3} - 6 x^{2} + 9 x = 0$

After:

$\displaystyle \left\{0, 3\right\}$

To get the solutions with polymonial equation with multiplicity, use roots()

In [33]:

eq = Eq(x**3 - 6*x**2 + 9*x,0)

func_comparator(eq, roots, x)

Before applying roots():

$\displaystyle x^{3} - 6 x^{2} + 9 x = 0$

After:

$\displaystyle \left\{ 0 : 1, \ 3 : 2\right\}$

Linear System

linsolve() is the solver to solve a linear system. A linear system can be defined using differet method.

list of equations

In [34]:

ls = [
    Eq(x + y + z - 1), 
    Eq(x + y + 2*z - 3)
]

func_comparator(ls,linsolve, x, z)

Before applying linsolve():

/usr/local/lib/python3.7/site-packages/sympy/core/relational.py:470: SymPyDeprecationWarning: 

Eq(expr) with rhs default to 0 has been deprecated since SymPy 1.5.
Use Eq(expr, 0) instead. See
https://github.com/sympy/sympy/issues/16587 for more info.

  deprecated_since_version="1.5"

$\displaystyle \left[ x + y + z - 1 = 0, \ x + y + 2 z - 3 = 0\right]$

After:

$\displaystyle \left\{\left( - y - 1, \ 2\right)\right\}$

List of equatons can be simplified into list of expressions if all of them are equal to zero。

In [35]:

ls = [
    x + y + z - 1, 
    x + y + 2*z - 3
]

Augmeted matrix form

In [36]:

M = Matrix(
    [
        [1, 1, 1, 1], 
        [1, 1, 2, 3]
    ]
)
    
func_comparator(M,linsolve, x, y, z)

Before applying linsolve():

$\displaystyle \left[\begin{matrix}1 & 1 & 1 & 1\\1 & 1 & 2 & 3\end{matrix}\right]$

After:

$\displaystyle \left\{\left( - y - 1, \ y, \ 2\right)\right\}$

$Ax = b$ form

In [37]:

M = Matrix(
    [
        [1, 1, 1, 1], 
        [1, 1, 2, 3]
    ]
)
ls = A, b = M[:, :-1], M[:, -1]

func_comparator(M,linsolve, x, y, z)

Before applying linsolve():

$\displaystyle \left[\begin{matrix}1 & 1 & 1 & 1\\1 & 1 & 2 & 3\end{matrix}\right]$

After:

$\displaystyle \left\{\left( - y - 1, \ y, \ 2\right)\right\}$

Non-linear Systems

solvset() is not capable of solving mutiple variable non-linear system. In this situation, use solve() instead.

In [38]:

eq_list = [
    Eq(x*y - 1, 0),
    Eq(x-2, 0)
]

func_comparator(eq_list, solve, x, y)

Before applying solve():

$\displaystyle \left[ x y - 1 = 0, \ x - 2 = 0\right]$

After:

$\displaystyle \left[ \left( 2, \ \frac{1}{2}\right)\right]$

Differential equition

First define needed function symbol variables.

In [39]:

f = symbols('f', cls = Function)

Then define the differential equation. For example, $f''(x) - 2f'(x) + f(x) = \sin(x)$

In [40]:

diffeq = Eq(f(x).diff(x, 2) - 2*f(x).diff(x) + f(x), sin(x))

Then use dsolve() to solve the equation.

In [41]:

func_comparator(diffeq, dsolve, f(x))

Before applying dsolve():

$\displaystyle f{\left(x \right)} - 2 \frac{d}{d x} f{\left(x \right)} + \frac{d^{2}}{d x^{2}} f{\left(x \right)} = \sin{\left(x \right)}$

After:

$\displaystyle f{\left(x \right)} = \left(C_{1} + C_{2} x\right) e^{x} + \frac{\cos{\left(x \right)}}{2}$

If the equation cannot be solved explicitly, dsolve() returns an implicit equation.

In [42]:

diffeq = Eq(f(x).diff(x)*(1 - sin(f(x))), 0)

func_comparator(diffeq, dsolve, f(x))

Before applying dsolve():

$\displaystyle \left(1 - \sin{\left(f{\left(x \right)} \right)}\right) \frac{d}{d x} f{\left(x \right)} = 0$

After:

$\displaystyle f{\left(x \right)} = C_{1}$

Reference

Sympy Documentation