[PATCH 1/2] Make Matrix.jacobian() work with any vectors
Ondrej Certik
calculating jacobians of row vectors. Now the whole algorithm
is implemented by a oneliner, everything else is just checks. A test was added
to make sure it works for any combination of row/column vectors.
Signed-off-by: Ondrej Certik <ond...@certik.cz>
---
sympy/matrices/matrices.py | 58 +++++++++++++++++++-------------
sympy/matrices/tests/test_matrices.py | 14 +++++++-
2 files changed, 47 insertions(+), 25 deletions(-)
diff --git a/sympy/matrices/matrices.py b/sympy/matrices/matrices.py
index 2e4e383..ad7c980 100644
--- a/sympy/matrices/matrices.py
+++ b/sympy/matrices/matrices.py
@@ -909,33 +909,43 @@ def cofactor(self, i, j, method="berkowitz"):
else:
return -1 * self.minorEntry(i, j, method)
- def jacobian(self, varlist):
+ def jacobian(self, X):
"""
Calculates the Jacobian matrix (derivative of a vectorial function).
- self is a vector of expression representing functions f_i(x_1, ...,
- x_n). varlist is the set of x_i's in order.
+ self ... a vector of expressions representing functions f_i(x_1, ...,
+ x_n).
+ X ...... is the set of x_i's in order, it can be a list or a Matrix
+
+ Both self and X can be a row or a column matrix in any order
+ (jacobian() should always work).
+
+ Example:
+ >>> from sympy import symbols, sin, cos
+ >>> rho, phi = symbols("rho phi")
+ >>> X = Matrix([rho*cos(phi), rho*sin(phi)])
+ >>> Y = Matrix([rho, phi])
+ >>> X.jacobian(Y)
+ [cos(phi), -rho*sin(phi)]
+ [sin(phi), rho*cos(phi)]
"""
- assert self.lines == 1
- m = self.cols
- if isinstance(varlist, Matrix):
- assert varlist.lines == 1
- n = varlist.cols
- elif isinstance(varlist, (list, tuple)):
- n = len(varlist)
- assert n > 0 # need to diff by something
- J = self.zeros((m, n)) # maintain subclass type
- for i in range(m):
- if isinstance(self[i], (float, int)):
- continue # constant function, jacobian row is zero
- try:
- tmp = self[i].diff(varlist[0]) # check differentiability
- J[i,0] = tmp
- except AttributeError:
- raise ValueError("Function %d is not differentiable" % i)
- for j in range(1,n):
- J[i,j] = self[i].diff(varlist[j])
- return J
+ if not isinstance(X, Matrix):
+ X = Matrix(X)
+ # Both X and self can be a row or a column matrix, so we need to make
+ # sure all valid combinations work, but everything else fails:
+ assert len(self.shape) == 2
+ assert len(X.shape) == 2
+ if self.shape[0] == 1:
+ n = self.shape[1]
+ else:
+ n = self.shape[0]
+ if X.shape[0] == 1:
+ assert X.shape[1] == n
+ else:
+ assert X.shape[0] == n
+
+ # n is the dimension of the matrix, computing the Jacobian is now easy:
+ return Matrix(n, n, lambda j, i: self[j].diff(X[i]))
def QRdecomposition(self):
"""
diff --git a/sympy/matrices/tests/test_matrices.py b/sympy/matrices/tests/test_matrices.py
index 98e93cf..6f71b52 100644
--- a/sympy/matrices/tests/test_matrices.py
+++ b/sympy/matrices/tests/test_matrices.py
@@ -1,6 +1,6 @@
from sympy import symbols, Matrix, eye, I, Symbol, Rational, wronskian, cos, \
sin, exp, hessian, sqrt, zeros, ones, randMatrix, Poly, S, pi, \
- integrate, oo, raises
+ integrate, oo, raises, trigsimp
from sympy.matrices.matrices import ShapeError, MatrixError
from sympy.printing import srepr
from sympy.utilities.pytest import XFAIL
@@ -957,3 +957,15 @@ def test_inv_iszerofunc():
A.col_swap(0,1)
for method in "GE", "LU":
assert A.inv(method, iszerofunc=lambda x: x==0) == A.inv("ADJ")
+
+def test_jacobian_metrics():
+ rho, phi = symbols("rho phi")
+ X = Matrix([rho*cos(phi), rho*sin(phi)])
+ Y = Matrix([rho, phi])
+ J = X.jacobian(Y)
+ assert J == X.jacobian(Y.T)
+ assert J == (X.T).jacobian(Y)
+ assert J == (X.T).jacobian(Y.T)
+ g = J.T*eye(J.shape[0])*J
+ g = g.applyfunc(trigsimp)
+ assert g == Matrix([[1, 0], [0, rho**2]])
--
1.6.2
Ondrej Certik
useful information about them.
---
examples/advanced/curvilinear_coordinates.py | 109 ++++++++++++++++++++++++++
examples/all.py | 1 +
2 files changed, 110 insertions(+), 0 deletions(-)
create mode 100644 examples/advanced/curvilinear_coordinates.py
diff --git a/examples/advanced/curvilinear_coordinates.py b/examples/advanced/curvilinear_coordinates.py
new file mode 100644
index 0000000..1436229
--- /dev/null
+++ b/examples/advanced/curvilinear_coordinates.py
@@ -0,0 +1,109 @@
+#!/usr/bin/env python
+
+"""
+This example shows how to work with coordinate transformations, curvilinear
+coordinates and a little bit with differential geometry.
+
+It takes polar, cylindrical, spherical, rotating disk coordinates and others
+and calculates all kinds of interesting properties, like Jacobian, metric
+tensor, Laplace operator, ...
+"""
+
+from sympy import var, sin, cos, pprint, Matrix, eye, trigsimp, Eq, \
+ Function, simplify, sinh, cosh
+
+def laplace(f, g_inv, g_det, X):
+ """
+ Calculates Laplace(f), using the inverse metric g_inv, the determinant of
+ the metric g_det, all in variables X.
+ """
+ r = 0
+ for i in range(len(X)):
+ for j in range(len(X)):
+ r += g_inv[i, j]*f.diff(X[i]).diff(X[j])
+ for sigma in range(len(X)):
+ for alpha in range(len(X)):
+ r += g_det.diff(X[sigma]) * g_inv[sigma, alpha] * \
+ f.diff(X[alpha]) / (2*g_det)
+ return r
+
+def transform(name, X, Y, g_correct=None):
+ """
+ Transforms from cartesian coordinates X to any curvilinear coordinates Y.
+
+ It printing useful information, like Jacobian, metric tensor, determinant
+ of metric, Laplace operator in the new coordinates, ...
+
+ g_correct ... if not None, it will be taken as the metric --- this is
+ useful if sympy's trigsimp() is not powerful enough to
+ simplify the metric so that it is usable for later
+ calculation. Leave it as None, only if the metric that
+ transform() prints is not simplified, you can help it by
+ specifying the correct one.
+ """
+ print "_"*80
+ print "Transformation:", name
+ for x, y in zip(X, Y):
+ pprint(Eq(y, x))
+ J = X.jacobian(Y)
+ print "Jacobian:"
+ pprint(J)
+ g = J.T*eye(J.shape[0])*J
+ g = g.applyfunc(trigsimp)
+ print "metric tensor g_{ij}:"
+ pprint(g)
+ if g_correct is not None:
+ g = g_correct
+ print "metric tensor g_{ij} specified by hand:"
+ pprint(g)
+ print "inverse metric tensor g^{ij}:"
+ g_inv = g.inv("ADJ")
+ g_inv = g_inv.applyfunc(simplify)
+ pprint(g_inv)
+ print "det g_{ij}:"
+ g_det = g.det()
+ pprint(g_det)
+ f = Function("f")(*list(Y))
+ print "Laplace:"
+ pprint(laplace(f, g_inv, g_det, Y))
+
+
+def main():
+ var("mu nu rho theta phi sigma tau a t x y z w")
+
+ transform("polar", Matrix([rho*cos(phi), rho*sin(phi)]), [rho, phi])
+
+ transform("cylindrical", Matrix([rho*cos(phi), rho*sin(phi), z]),
+ [rho, phi, z])
+
+ transform("spherical",
+ Matrix([rho*sin(theta)*cos(phi), rho*sin(theta)*sin(phi),
+ rho*cos(theta)]),
+ [rho, theta, phi],
+ # spherical metric needs a little help:
+ g_correct=Matrix([[1, 0, 0], [0, rho**2, 0],
+ [0, 0, rho**2*sin(theta)**2]])
+ )
+
+ transform("rotating disk",
+ Matrix([t, x*cos(w*t)-y*sin(w*t), x*sin(w*t)+y*cos(w*t), z]),
+ [t, x, y, z])
+
+ transform("parabolic",
+ Matrix([sigma*tau, (tau**2-sigma**2)/2]),
+ [sigma, tau])
+
+ # too complex:
+ #transform("bipolar",
+ # Matrix([a*sinh(tau)/(cosh(tau)-cos(sigma)),
+ # a*sin(sigma)/(cosh(tau)-cos(sigma))]),
+ # [sigma, tau]
+ # )
+
+ transform("elliptic",
+ Matrix([a*cosh(mu)*cos(nu), a*sinh(mu)*sin(nu)]),
+ [mu, nu]
+ )
+
+if __name__ == "__main__":
+ main()
diff --git a/examples/all.py b/examples/all.py
index bacfc7e..46e460b 100755
--- a/examples/all.py
+++ b/examples/all.py
@@ -62,6 +62,7 @@
"advanced.pidigits",
"advanced.qft",
"advanced.relativity",
+ "advanced.curvilinear_coordinates",
]
WINDOWED_EXAMPLES = [
--
1.6.2
Fabian Seoane
Ondrej Certik
O.