Juliawe have
Juliasave for ** becoming ^ this is the same
diff(cos(x), x)
diff(exp(x^2), x)
Julia:@vars m n a b
sin(x)
Julia:integrate(cos(x), x)
Note that SymPy does not include the constant of integration. If you want it, you can add one yourself, or rephrase your problem as a differential equation and use dsolve to solve it, which does add the constant (see :ref:tutorial-dsolve).
Quick Tip:
@@ -266,11 +272,11 @@In Julia:
expr = integrate(x^x, x)In
Julia:the
Integralclass is not exported, so it must be qualified:expr = sympy.Integral(log(x)^2, x) expr-\begin{equation*}\int \log{\left (x \right )}^{2}\, dx\end{equation*}+\begin{equation*}\int \log{\left(x \right)}^{2}\, dx\end{equation*}expr.doit()-\begin{equation*}x \log{\left (x \right )}^{2} - 2 x \log{\left (x \right )} + 2 x\end{equation*}+\begin{equation*}x \log{\left(x \right)}^{2} - 2 x \log{\left(x \right)} + 2 x\end{equation*}
integrateuses powerful algorithms that are always improving to compute both definite and indefinite integrals, including heuristic pattern matching type algorithms, a partial implementation of theRisch algorithm <http://en.wikipedia.org/wiki/Risch_algorithm>, and an algorithm usingMeijer G-functions <http://en.wikipedia.org/wiki/Meijer_g-function>that is useful for computing integrals in terms of special functions, especially definite integrals. Here is a sampling of some of the power ofintegrate.>>> integ = Integral((x**4 + x**2*exp(x) - x**2 - 2*x*exp(x) - 2*x - ... exp(x))*exp(x)/((x - 1)**2*(x + 1)**2*(exp(x) + 1)), x) @@ -327,7 +333,7 @@In
Julia:integ = sympy.Integinteg.doit()-\begin{equation*}\log{\left (e^{x} + 1 \right )} + \frac{e^{x}}{x^{2} - 1}\end{equation*}+\begin{equation*}\log{\left(e^{x} + 1 \right)} + \frac{e^{x}}{x^{2} - 1}\end{equation*}integ = sympy.Integral(sin(x^2), x) integ.doit()@@ -337,7 +343,7 @@In
Julia:integ = sympy.Integinteg = sympy.Integral(x^y*exp(-x), (x, 0, oo)) integ.doit()-\begin{equation*}\begin{cases} \Gamma\left(y + 1\right) & \text{for}\: - \Re{\left(y\right)} < 1 \\\int\limits_{0}^{\infty} x^{y} e^{- x}\, dx & \text{otherwise} \end{cases}\end{equation*}+\begin{equation*}\begin{cases} \Gamma\left(y + 1\right) & \text{for}\: \operatorname{re}{\left(y\right)} > -1 \\\int\limits_{0}^{\infty} x^{y} e^{- x}\, dx & \text{otherwise} \end{cases}\end{equation*}This last example returned a
Piecewiseexpression because the integral does not converge unless $\Re(y) > 1.$Limits
SymPy can compute symbolic limits with the
$$~ \lim_{x\to x_0} f(x) @@ -360,7 +366,7 @@limitfunction. The syntax to computeIn
Julia:limit(sin(x)/x, x,In
Julia:expr = x^2/exp(x) expr.subs(x, oo)-\begin{equation*}\mathrm{NaN}\end{equation*}+\begin{equation*}\text{NaN}\end{equation*}limit(expr, x, oo)@@ -376,7 +382,7 @@In
Julia:expr = x^2/exp(x)In
Julia:expr = sympy.Limit((cos(x) - 1)/x, x, 0) expr-\begin{equation*}\lim_{x \to 0^+}\left(\frac{\cos{\left (x \right )} - 1}{x}\right)\end{equation*}+\begin{equation*}\lim_{x \to 0^+}\left(\frac{\cos{\left(x \right)} - 1}{x}\right)\end{equation*}expr.doit()@@ -449,13 +455,13 @@In
Julia:
f, g = symbols("f g", cls=sympy.Function) sympy.differentiate_finite(f(x)*g(x))-\begin{equation*}- f{\left (x - \frac{1}{2} \right )} g{\left (x - \frac{1}{2} \right )} + f{\left (x + \frac{1}{2} \right )} g{\left (x + \frac{1}{2} \right )}\end{equation*}+\begin{equation*}- f{\left(x - \frac{1}{2} \right)} g{\left(x - \frac{1}{2} \right)} + f{\left(x + \frac{1}{2} \right)} g{\left(x + \frac{1}{2} \right)}\end{equation*}(The functions
fandgcan also be created with the command@symfuns f g, using the@symfunsmacro.)If we want to expand the intermediate derivative we may pass the flag
evaluate=True:>>> differentiate_finite(f(x)*g(x), evaluate=True) (-f(x - 1/2) + f(x + 1/2))⋅g(x) + (-g(x - 1/2) + g(x + 1/2))⋅f(x)In
Julia:sympy.differentiate_finite(f(x)*g(x), evaluate=true)-\begin{equation*}\left(- f{\left (x - \frac{1}{2} \right )} + f{\left (x + \frac{1}{2} \right )}\right) g{\left (x \right )} + \left(- g{\left (x - \frac{1}{2} \right )} + g{\left (x + \frac{1}{2} \right )}\right) f{\left (x \right )}\end{equation*}+\begin{equation*}\left(- f{\left(x - \frac{1}{2} \right)} + f{\left(x + \frac{1}{2} \right)}\right) g{\left(x \right)} + \left(- g{\left(x - \frac{1}{2} \right)} + g{\left(x + \frac{1}{2} \right)}\right) f{\left(x \right)}\end{equation*}This form however does not respect the product rule.
If you already have a
Derivativeinstance, you can use theas_finite_differencemethod to generate approximations of the derivative to arbitrary order:>>> f = Function('f') >>> dfdx = f(x).diff(x) @@ -465,7 +471,7 @@In
Julia:f = sympy.Function( dfdx = f(x).diff(x) dfdx.as_finite_difference()-\begin{equation*}- f{\left (x - \frac{1}{2} \right )} + f{\left (x + \frac{1}{2} \right )}\end{equation*}+\begin{equation*}- f{\left(x - \frac{1}{2} \right)} + f{\left(x + \frac{1}{2} \right)}\end{equation*}here the first order derivative was approximated around x using a minimum number of points (2 for 1st order derivative) evaluated equidistantly using a step-size of 1. We can use arbitrary steps (possibly containing symbolic expressions):
>>> f = Function('f') >>> d2fdx2 = f(x).diff(x, 2) @@ -480,7 +486,7 @@In
Julia:f = sympy.Function( h = sympy.Symbol("h") d2fdx2.as_finite_difference([-3*h,-h,2*h])-\begin{equation*}\frac{f{\left (- 3 h \right )}}{5 h^{2}} - \frac{f{\left (- h \right )}}{3 h^{2}} + \frac{2 f{\left (2 h \right )}}{15 h^{2}}\end{equation*}+\begin{equation*}\frac{f{\left(- 3 h \right)}}{5 h^{2}} - \frac{f{\left(- h \right)}}{3 h^{2}} + \frac{2 f{\left(2 h \right)}}{15 h^{2}}\end{equation*}If you are just interested in evaluating the weights, you can do so manually:
>>> finite_diff_weights(2, [-3, -1, 2], 0)[-1][-1] [1/5, -1/3, 2/15]diff --git a/examples/matrices.html b/examples/matrices.html index 5a8bfc35..1b677c3a 100644 --- a/examples/matrices.html +++ b/examples/matrices.html @@ -133,17 +133,8 @@In
Julia:
In
Julia, theMatrixconstructor is not exported, so must be qualified. Here we avoid thte vector of row vectors above:sympy.Matrix([1 -1; 3 4; 0 2])--\[\left[ \begin{array}{rr}1&-1\\3&4\\0&2\end{array}\right]\]-However, through the magic of
PyCall, such matrices are converted intoJuliamatrices, of typeArray{Sym}, so the familiar matrix operations forJuliausers are available.In fact, the above could be done in the more
Julian manner throughSym[1 -1; 3 4; 0 2]--\[\left[ \begin{array}{rr}1&-1\\3&4\\0&2\end{array}\right]\]-using an annotation to ensure the type. Alternatively, through promotion, just a single symbolic object will result in the same:
[Sym(1) -1; 3 4; 0 2]--\[\left[ \begin{array}{rr}1&-1\\3&4\\0&2\end{array}\right]\]-To make it easy to make column vectors, a list of elements is considered to be a column vector.
>>> Matrix([1, 2, 3]) ⎡1⎤ ⎢ ⎥ @@ -151,16 +142,10 @@In
Julia:⎢ ⎥ ⎣3⎦
In
Julia:For ths use,
sympy.Matrixdoes work, but again its usage is discouraged:sympy.Matrix([1, 2, 3])--\[\left[ \begin{array}{r}1\\2\\3\end{array}\right]\]-
Again, this is converted into a
Vector{Sym}object or entered directly:Sym[1,2,3]--\[ \left[ \begin{array}{r}1\\2\\3\end{array} \right] \]-And again:
Using
sympy.Matrixis strongly discouraged.Matrices are manipulated just like any other object in SymPy or Python.
>>> M = Matrix([[1, 2, 3], [3, 2, 1]]) >>> N = Matrix([0, 1, 1]) @@ -174,11 +159,6 @@In
Julia:
M = Sym[1 2 3; 3 2 1] N = Sym[0, 1, 1] M*N--2-element Array{Any,1}: - 5 - 3-One important thing to note about SymPy matrices is that, unlike every other object in SymPy, they are mutable. This means that they can be modified in place, as we will see below. The downside to this is that
Matrixcannot be used in places that require immutability, such as inside other SymPy expressions or as keys to dictionaries. If you need an immutable version ofMatrix, useImmutableMatrix.In
Julia:A distinction is made between
ImmutableMatrixand a mutable one. Mutable ones are mapped toJuliaarrays, immutable ones are left as a symbolic object of typeSymMatrix. The usual infix mathematical operations (but not dot broadcasting), 0-based indexing, and dot call syntax for methods are used with these objects.Basic Operations
Shape
Here are some basic operations on
Matrix. To get the shape of a matrix useshape>>> M = Matrix([[1, 2, 3], [-2, 0, 4]]) >>> M ⎡1 2 3⎤ @@ -188,9 +168,6 @@In
Julia:(2, 3)
In
Julia:M = Sym[1 2 3; -2 0 4] M--\[\left[ \begin{array}{rrr}1&2&3\\-2&0&4\end{array}\right]\]-M.shape(2, 3)@@ -210,12 +187,9 @@In
Julia:
M.row(0) M.col(-1)--\[\left[ \begin{array}{r}3\\4\end{array}\right]\]-The more familiar counterparts would be:
M[1,:], M[:, end]-(Sym[1, 2, 3], Sym[3, 4])+(SymPy.Sym[SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)), SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000146026dd0)), SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000146405b30))], SymPy.Sym[SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000146405b30)), SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000146d01530))])Deleting and Inserting Rows and Columns
To delete a row or column, use
row_delorcol_del. These operations will modify the Matrix in place.>>> M.col_del(0) >>> M @@ -227,18 +201,9 @@In
Julia:[2 3]
In
Julia:These methods do not work on
Array{Sym}objects, useJulia'sindexing notation to remove a row or column.However, these methods do work on the
ImmutableMatrixclass:M = sympy.ImmutableMatrix([1 2 3; -2 0 4]) # avoid vector of row vector construction M.col_del(0)--\begin{equation*}\left[\begin{matrix}2 & 3\\0 & 4\end{matrix}\right]\end{equation*}-M.row_del(1)--\begin{equation*}\left[\begin{matrix}1 & 2 & 3\end{matrix}\right]\end{equation*}-!!! Alert
For older versions of sympy, the following did not work (using symbolic values as matrix entries without reverting to their PyObjects had shape issues) but this should work now:
@vars x sympy.ImmutableMatrix([x 1; 1 x])--\begin{equation*}\left[\begin{matrix}x & 1\\1 & x\end{matrix}\right]\end{equation*}-TODO
This is a mess. See issue 6992.
To insert rows or columns, use
row_insertorcol_insert. These operations do not operate in place.>>> M [2 3] @@ -255,14 +220,8 @@In
Julia:These methods do not work on
In
Julia:M M = M.row_insert(1, Sym[0 4]) M--\begin{equation*}\left[\begin{matrix}1 & 2 & 3\\-2 & 0 & 4\end{matrix}\right]\end{equation*}-M = M.col_insert(0, Matrix([1, -2])) M--\begin{equation*}\left[\begin{matrix}1 & 2 & 3\\-2 & 0 & 4\end{matrix}\right]\end{equation*}-Unless explicitly stated, the methods mentioned below do not operate in place. In general, a method that does not operate in place will return a new
Matrixand a method that does operate in place will returnNone.In
JuliaThis would be the case for the immutable matrices.
Basic Methods
As noted above, simple operations like addition and multiplication are done just by using
+,*, and**. To find the inverse of a matrix, just raise it to the-1power.>>> M = Matrix([[1, 3], [-2, 3]]) >>> N = Matrix([[0, 3], [0, 7]]) >>> M + N @@ -292,75 +251,25 @@In
Julia:MIn
Julia:M = Sym[1 3; -2 3] M1 = Sym[0 3; 0 7] M + M1--\[\left[ \begin{array}{rr}1&6\\-2&10\end{array}\right]\]-M*M1--2×2 Array{Any,2}: - 0 24 - 0 15-3*M--\[\left[ \begin{array}{rr}3&9\\-6&9\end{array}\right]\]-M^2--2×2 Array{Any,2}: - -5 12 - -8 3-M^-1--\[\left[ \begin{array}{rr}\frac{1}{3}&- \frac{1}{3}\\\frac{2}{9}&\frac{1}{9}\end{array}\right]\]-M1^-1-PyError ($(Expr(:escape, :(ccall(#= /Users/verzani/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:44 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw)))))-NonInvertibleMatrixError('Matrix det == 0; not invertible.') - File "/Users/verzani/.julia/conda/3/lib/python3.7/site-packages/sympy/matrices/matrices.py", line 3672, in inv - return self._eval_inverse(**kwargs) - File "/Users/verzani/.julia/conda/3/lib/python3.7/site-packages/sympy/matrices/dense.py", line 265, in _eval_inverse - rv = M.inverse_GE(iszerofunc=iszerofunc) - File "/Users/verzani/.julia/conda/3/lib/python3.7/site-packages/sympy/matrices/matrices.py", line 3587, in inverse_GE - raise NonInvertibleMatrixError("Matrix det == 0; not invertible.") - +PyCall.PyError("\$(Expr(:escape, :(ccall(#= /Users/verzani/.julia/packages/PyCall/zqDXB/src/pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))", PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x00007fe06ca00950), PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x000000015363ec90), PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000153678280)) The above (except for the inverses) are using generic
Juliadefinitions. For immutable matrices, we would have:M = sympy.ImmutableMatrix([1 3; -2 3]) M1 = sympy.ImmutableMatrix([0 3; 0 7]) M + M1--\begin{equation*}\left[\begin{matrix}1 & 6\\-2 & 10\end{matrix}\right]\end{equation*}-M*M1--\begin{equation*}\left[\begin{matrix}0 & 24\\0 & 15\end{matrix}\right]\end{equation*}-3*M--\begin{equation*}\left[\begin{matrix}3 & 9\\-6 & 9\end{matrix}\right]\end{equation*}-M^2--\begin{equation*}\left(\left[\begin{matrix}1 & 3\\-2 & 3\end{matrix}\right]\right)^{2}\end{equation*}-M^-1--\begin{equation*}\left[\begin{matrix}\frac{1}{3} & - \frac{1}{3}\\\frac{2}{9} & \frac{1}{9}\end{matrix}\right]\end{equation*}-M1^-1-PyError ($(Expr(:escape, :(ccall(#= /Users/verzani/.julia/packages/PyCall/ttONZ/src/pyfncall.jl:44 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw)))))-NonInvertibleMatrixError('Matrix det == 0; not invertible.') - File "/Users/verzani/.julia/conda/3/lib/python3.7/site-packages/sympy/matrices/matrices.py", line 3672, in inv - return self._eval_inverse(**kwargs) - File "/Users/verzani/.julia/conda/3/lib/python3.7/site-packages/sympy/matrices/dense.py", line 265, in _eval_inverse - rv = M.inverse_GE(iszerofunc=iszerofunc) - File "/Users/verzani/.julia/conda/3/lib/python3.7/site-packages/sympy/matrices/matrices.py", line 3587, in inverse_GE - raise NonInvertibleMatrixError("Matrix det == 0; not invertible.") - +PyCall.PyError("\$(Expr(:escape, :(ccall(#= /Users/verzani/.julia/packages/PyCall/zqDXB/src/pyfncall.jl:43 =# @pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, pyargsptr, kw))))", PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x00007fe06ca00950), PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000153c0a7c0), PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x00000001536819b0)) @@ -379,13 +288,7 @@
In
Julia:M = Sym[1 3; -2 3] ⎣3 6⎦In
Julia:M = Sym[1 2 3; 4 5 6] M--\[\left[ \begin{array}{rrr}1&2&3\\4&5&6\end{array}\right]\]-M.T--\[\left[ \begin{array}{rr}1&4\\2&5\\3&6\end{array}\right]\]-Matrix Constructors
Several constructors exist for creating common matrices. To create an identity matrix, use
eye. The commandeye(n)will create ann x nidentity matrix:>>> eye(3) ⎡1 0 0⎤ ⎢ ⎥ @@ -405,9 +308,6 @@In
Julia:
sympy.eye(3) sympy.eye(4)--\[\left[ \begin{array}{rrrr}1&0&0&0\\0&1&0&0\\0&0&1&0\\0&0&0&1\end{array}\right]\]-To create a matrix of all zeros, use
zeros.zeros(n, m)creates ann x mmatrix of0s.>>> zeros(2, 3) ⎡0 0 0⎤ ⎢ ⎥ @@ -416,17 +316,8 @@In
Julia:
zeros is exported but the method expects a symbolic first argument. Either qualify it:
sympy.zeros(2, 3)--\[\left[ \begin{array}{rrr}0&0&0\\0&0&0\end{array}\right]\]-or create a symbolic first value:
zeros(Sym(2), 3)--\[\left[ \begin{array}{rrr}0&0&0\\0&0&0\end{array}\right]\]-or use the
Juliaconstructor:zeros(Sym, 2, 3)--\[\left[ \begin{array}{rrr}0&0&0\\0&0&0\end{array}\right]\]-Similarly,
onescreates a matrix of ones.>>> ones(3, 2) ⎡1 1⎤ ⎢ ⎥ @@ -437,9 +328,6 @@In
Julia:
Similarly with
ones:sympy.ones(3, 2)--\[\left[ \begin{array}{rr}1&1\\1&1\\1&1\end{array}\right]\]-To create diagonal matrices, use
diag. The arguments todiagcan be either numbers or matrices. A number is interpreted as a1 x 1matrix. The matrices are stacked diagonally. The remaining elements are filled with0\ s.>>> diag(1, 2, 3) ⎡1 0 0⎤ ⎢ ⎥ @@ -462,21 +350,12 @@In
Julia:
similarly with
diag:sympy.diag(1, 2, 3)--\[\left[ \begin{array}{rrr}1&0&0\\0&2&0\\0&0&3\end{array}\right]\]-sympy.diag(-1, sympy.ones(2, 2), sympy.Matrix([5, 7, 5]))--\[\left[ \begin{array}{rrrr}-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{array}\right]\]-
The first one, could also use
Julia'sdiagmfunction:using LinearAlgebra diagm(0 => Sym[1,2,3])--\[\left[ \begin{array}{rrr}1&0&0\\0&2&0\\0&0&3\end{array}\right]\]-Advanced Methods
Determinant
To compute the determinant of a matrix, use
det.>>> M = Matrix([[1, 0, 1], [2, -1, 3], [4, 3, 2]]) >>> M ⎡1 0 1⎤ @@ -489,26 +368,11 @@In
Julia:
In
Julia:using LinearAlgebra M = Sym[1 0 1; 2 -1 3; 4 3 2] M--\[\left[ \begin{array}{rrr}1&0&1\\2&-1&3\\4&3&2\end{array}\right]\]-M.det()--\begin{equation*}-1\end{equation*}-Let
@vars x A = Sym[x 1; 1 x]--\[\left[ \begin{array}{rr}x&1\\1&x\end{array}\right]\]-Then we can compute the determinant using
Julia's generic implementation:det(A)--\begin{equation*}x^{2} - 1\end{equation*}-or using SymPy's:
A.det()--\begin{equation*}x^{2} - 1\end{equation*}-The answer is identical, though not necessarily being done in a similar manner.
RREF
To put a matrix into reduced row echelon form, use
rref.rrefreturns a tuple of two elements. The first is the reduced row echelon form, and the second is a tuple of indices of the pivot columns.>>> M = Matrix([[1, 0, 1, 3], [2, 3, 4, 7], [-1, -3, -3, -4]]) >>> M ⎡1 0 1 3 ⎤ @@ -524,12 +388,9 @@In
Julia:using LinearAlgebra ⎝⎣0 0 0 0 ⎦ ⎠In
Julia:M = Sym[1 0 1 3; 2 3 4 7; -1 -3 -3 -4] M--\[\left[ \begin{array}{rrrr}1&0&1&3\\2&3&4&7\\-1&-3&-3&-4\end{array}\right]\]-M.rref()-(Sym[Sym(PyObject 1) Sym(PyObject 0) Sym(PyObject 1) Sym(PyObject 3); Sym(PyObject 0) Sym(PyObject 1) Sym(PyObject 2/3) Sym(PyObject 1/3); Sym(PyObject 0) Sym(PyObject 0) Sym(PyObject 0) Sym(PyObject 0)], (0, 1))+(SymPy.Sym[SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)) SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30)) SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)) SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000146405b30)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30)) SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)) SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x000000015367fa10)) SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000153699b30)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30)) SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30)) SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30)) SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30))], (0, 1))Note
The first element of the tuple returned by
rrefis of typeMatrix. The second is of typetuple.Nullspace
To find the nullspace of a matrix, use
nullspace.nullspacereturns alistof column vectors that span the nullspace of the matrix.>>> M = Matrix([[1, 2, 3, 0, 0], [4, 10, 0, 0, 1]]) @@ -552,15 +413,12 @@In
Julia:
M = Sym[1 2 3 0 0; 4 10 0 0 1] M--\[\left[ \begin{array}{rrrrr}1&2&3&0&0\\4&10&0&0&1\end{array}\right]\]-M.nullspace()-3-element Array{Array{Sym,2},1}: - [-15; 6; 1; 0; 0] - [0; 0; 0; 1; 0] - [1; -1/2; 0; 0; 1]+3-element Array{Array{SymPy.Sym,2},1}: + [SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000153c1ebf0)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000146405470)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30))] + [SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30))] + [SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000153c1ecb0)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650))]Columnspace
To find the columnspace of a matrix, use
columnspace.columnspacereturns alistof column vectors that span the columnspace of the matrix.>>> M = Matrix([[1, 1, 2], [2 ,1 , 3], [3 , 1, 4]]) >>> M @@ -580,14 +438,11 @@In
Julia:
M = Sym[1 1 2; 2 1 3; 3 1 4] M--\[\left[ \begin{array}{rrr}1&1&2\\2&1&3\\3&1&4\end{array}\right]\]-M.columnspace()-2-element Array{Array{Sym,2},1}: - [1; 2; 3] - [1; 1; 1]+2-element Array{Array{SymPy.Sym,2},1}: + [SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000146026dd0)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000146405b30))] + [SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650))]Eigenvalues, Eigenvectors, and Diagonalization
To find the eigenvalues of a matrix, use
eigenvals.eigenvalsreturns a dictionary ofeigenvalue:algebraic multiplicitypairs (similar to the output of :ref:roots <tutorial-roots>).>>> M = Matrix([[3, -2, 4, -2], [5, 3, -3, -2], [5, -2, 2, -2], [5, -2, -3, 3]]) >>> M @@ -602,15 +457,12 @@In
Julia:{-2: 1, 3: 1, 5: 2}
In
Julia:M = Sym[3 -2 4 -2; 5 3 -3 -2; 5 -2 2 -2; 5 -2 -3 3] M--\[\left[ \begin{array}{rrrr}3&-2&4&-2\\5&3&-3&-2\\5&-2&2&-2\\5&-2&-3&3\end{array}\right]\]-M.eigenvals()Dict{Any,Any} with 3 entries: - 3 => 1 - -2 => 1 - 5 => 2+ SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000146405b30)) => 1 + SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x00000001222c57d0)) => 1 + SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x000000014602fb90)) => 2This means that
Mhas eigenvalues -2, 3, and 5, and that the eigenvalues -2 and 3 have algebraic multiplicity 1 and that the eigenvalue 5 has algebraic multiplicity 2.To find the eigenvectors of a matrix, use
eigenvects.eigenvectsreturns a list of tuples of the form(eigenvalue:algebraic multiplicity, [eigenvectors]).>>> M.eigenvects() ⎡⎛ ⎡⎡0⎤⎤⎞ ⎛ ⎡⎡1⎤⎤⎞ ⎛ ⎡⎡1⎤ ⎡0 ⎤⎤⎞⎤ @@ -627,15 +479,12 @@In
Julia:
M.eigenvects()-3-element Array{Tuple{Sym,Int64,Array{Array{Sym,2},1}},1}: - (-2, 1, [[0; 1; 1; 1]]) - (3, 1, [[1; 1; 1; 1]]) - (5, 2, [[1; 1; 1; 0], [0; -1; 0; 1]])+3-element Array{Tuple{SymPy.Sym,Int64,Array{Array{SymPy.Sym,2},1}},1}: + (SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x00000001222c57d0)), 1, [[SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650))]]) + (SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000146405b30)), 1, [[SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650))]]) + (SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x000000014602fb90)), 2, [[SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30))], [SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x00000001464055f0)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145c93b30)); SymPy.Sym(PyCall.PyObject(Ptr{PyCall.PyObject_struct} @0x0000000145be0650))]])compare with
eigvecs(M)--\[\left[ \begin{array}{rrrr}0&1&1&0\\1&1&1&-1\\1&1&1&0\\1&1&0&1\end{array}\right]\]-This shows us that, for example, the eigenvalue 5 also has geometric multiplicity 2, because it has two eigenvectors. Because the algebraic and geometric multiplicities are the same for all the eigenvalues,
Mis diagonalizable.To diagonalize a matrix, use
diagonalize.diagonalizereturns a tuple(P, D), whereDis diagonal andM = PDP^{-1}.>>> P, D = M.diagonalize() >>> P ⎡0 1 1 0 ⎤ @@ -665,21 +514,8 @@In
Julia:True
In
Julia:P, D = M.diagonalize() P--\[\left[ \begin{array}{rrrr}0&1&1&0\\1&1&1&-1\\1&1&1&0\\1&1&0&1\end{array}\right]\]-D--\[\left[ \begin{array}{rrrr}-2&0&0&0\\0&3&0&0\\0&0&5&0\\0&0&0&5\end{array}\right]\]-P*D*P^-1--4×4 Array{Any,2}: - 3 -2 4 -2 - 5 3 -3 -2 - 5 -2 2 -2 - 5 -2 -3 3-P*D*P^-1 == Mtrue@@ -695,9 +531,6 @@In
Julia:
lambda = symbols("lambda") p = M.charpoly(lambda) factor(p)--\begin{equation*}\left(\lambda - 5\right)^{2} \left(\lambda - 3\right) \left(\lambda + 2\right)\end{equation*}-TODO
Add an example for
jordan_form, once it is fully implemented.Possible Issues
Zero Testing
If your matrix operations are failing or returning wrong answers, the common reasons would likely be from zero testing. If there is an expression not properly zero-tested, it can possibly bring issues in finding pivots for gaussian elimination, or deciding whether the matrix is inversible, or any high level functions which relies on the prior procedures.
Currently, the SymPy's default method of zero testing
_iszerois only guaranteed to be accurate in some limited domain of numerics and symbols, and any complicated expressions beyond its decidability are treated asNone, which behaves similarly to logicalFalse.The list of methods using zero testing procedures are as followings.
echelon_form,is_echelon,rank,rref,nullspace,eigenvects,inverse_ADJ,inverse_GE,inverse_LU,LUdecomposition,LUdecomposition_Simple,LUsolveThey have property
iszerofuncopened up for user to specify zero testing method, which can accept any function with single input and boolean output, while being defaulted with_iszero.Here is an example of solving an issue caused by undertested zero. [#zerotestexampleidea-fn]_ [#zerotestexamplediscovery-fn]_
>>> from sympy import * >>> q = Symbol("q", positive = True) diff --git a/src/SymPy.jl b/src/SymPy.jl index 8e24f718..a4bec7d3 100644 --- a/src/SymPy.jl +++ b/src/SymPy.jl @@ -49,7 +49,7 @@ export SymMatrix, SymFunction export PI, IM, oo, zoo, True, False export N, subs -export sympy, import_from#, import_sympy +export sympy, sympy_core, sympy_matrices, import_from#, import_sympy export free_symbols @@ -72,8 +72,23 @@ include("plot_recipes.jl") ################################################## pynull() = PyCall.PyNULL() +""" + sympy + +variable from `pyimport("sympy")` +""" const sympy = PyCall.PyNULL() +""" + sympy_core + +variable from `pyimport("sympy.core")` +""" const sympy_core = PyCall.PyNULL() +""" + sympy_matrices + +variable from `pyimport("sympy.matrices")` +""" const sympy_matrices = PyCall.PyNULL() const mpmath = PyCall.PyNULL() const combinatorics = PyCall.PyNULL() diff --git a/src/mathops.jl b/src/mathops.jl index 5b8e3bd6..e0d8a3ca 100644 --- a/src/mathops.jl +++ b/src/mathops.jl @@ -23,12 +23,12 @@ inv(x::Sym) = x.__pow__(Sym(-1)) # special case Boolean; issue 351 # promotion for Boolean here is to 0 or 1, not False, True +(x::Bool, y::SymbolicObject) = Sym(Int(x)).__add__(y) -+(x::SymbolicObject, y::Bool) = x.__add__(Int(y)) *(x::Bool, y::SymbolicObject) = Sym(Int(x)).__mul__(y) -*(x::SymbolicObject, y::Bool) = x.__mul__(Int(y)) -(x::Bool, y::SymbolicObject) = Sym(Int(x)).__sub__(y) --(x::SymbolicObject, y::Bool) = x.__sub__(Int(y)) /(x::Bool, y::SymbolicObject) = Sym(Int(x)).__div__(y) -/(x::SymbolicObject, y::Bool) = x.__div__(Int(y)) ^(x::Bool, y::SymbolicObject) = Sym(Int(x)).__pow__(y) ++(x::SymbolicObject, y::Bool) = x.__add__(Int(y)) +*(x::SymbolicObject, y::Bool) = x.__mul__(Int(y)) +-(x::SymbolicObject, y::Bool) = x.__sub__(Int(y)) +/(x::SymbolicObject, y::Bool) = x.__div__(Int(y)) ^(x::SymbolicObject, y::Bool) = x.__pow__(Int(y)) diff --git a/test/tests.jl b/test/tests.jl index 71649686..170a99ae 100644 --- a/test/tests.jl +++ b/test/tests.jl @@ -655,10 +655,10 @@ end ## Issue 351 booleans and arithmetic operations - @test 1 + true == 2 == true + 1 - @test 1 - true == 0 == true - 1 - @test 1 * true == 1 == true * 1 - @test 1/true == 1 == true/1 - @test true^1 == 1 == 1^true + @test Sym(1) + true == Sym(2) == true + Sym(1) + @test Sym(1) - true == Sym(0) == true - Sym(1) + @test Sym(1) * true == Sym(1) == true * Sym(1) + @test Sym(1) / true == Sym(1) == true / Sym(1) + @test true^Sym(1) == Sym(1) == Sym(1)^true end