A computer algebra system written in pure Python http://sympy.org/ . To get started to with contributing https://github.com/sympy/sympy/wiki/Introduction-to-contributing
oscarbenjamin on master
core/containers: raising error … create tuple of comprehension s… recursively convert list argume… and 18 more (compare)
Hi guys, question about root finding with symbolic coefficients. Is it possible to get symbolic results for something like the following:
import sympy as sp
t = sp.Symbol('t')
P = sp.Poly(t**2 + 2*t + sp.log(2))
P.all_roots()
This does not work, as I get sympy.polys.polyerrors.PolynomialError: only univariate polynomials are allowed
P = sp.Poly(t**2 + 2*t + sp.log(2), t, domain='RR')
works for numerical values. Looking at the docs, I thought something like P = sp.Poly(t**2 + 2*t + sp.log(2), t, domain='QQ<log(2)>')
might work, but this triggers a sympy.polys.polyerrors.NotAlgebraic: log(2) doesn't seem to be an algebraic element
QQ<log(2)>
does not work because log(2)
is not algebraic. One must use something like QQ(log(2))
instead. Or, simply give t
as the sole generator; then log(2)
will automatically become a constant.In [26]: p = Poly(t**2 + 2*t + log(2), t)
In [27]: roots(p)
Out[27]:
⎧ ____________ ____________ ⎫
⎨-1 - ╲╱ 1 - log(2) : 1, -1 + ╲╱ 1 - log(2) : 1⎬
⎩ ⎭
Sum
, replicating in a "single command" things likeGeneralSumOfSquares
class used here. Is there a paper I can reference/explanation, etc?
I'm having trouble learning the diffgeom
module. I'm trying to compute covariant derivatives in spherical coordinates and I feel like I'm missing something.
Here's what I'm working with:
from sympy import *
from sympy.abc import x, y, z, r, theta, phi
from sympy.diffgeom import CoordSystem, Manifold, Patch, TensorProduct as TP
R3 = Manifold("R3", 3)
S = Patch("S", R3)
relations = {
("Car3D", "Sph") : Lambda(
(x, y, z),
Matrix([
sqrt(x**2 + y**2 + z**2),
acos(z / sqrt(x**2 + y**2 + z**2)),
atan2(y, x),
])
),
("Sph", "Car3D"): Lambda(
(r, theta, phi),
Matrix([r * sin(theta) * cos(phi), r * sin(theta) * sin(phi), r * cos(theta)])
)
}
Car3D = CoordSystem("Car3D", S, [x, y, z], relations)
Sph = CoordSystem("Sph", S, [r, theta, phi], relations)
I can compute the metric tensor in spherical coordinates by doing the following:
g = (Car3D.jacobian(Sph) @ Car3D.jacobian(Sph).T).simplify()
However, this computes the metric tensor with respect to $x, y, z$, while I'd like to write the metric tensor with respect to $r, \theta, \phi$:
$\left[\begin{matrix}1 & 0 & 0\\0 & \frac{1}{x^{2} + y^{2} + z^{2}} & 0\\0 & 0 & \frac{1}{x^{2} + y^{2}}\end{matrix}\right]$
I can do a substitution, but it's kind of awkward:
x, y, z = Car3D.symbols
r, theta, phi = Sph.symbols
forward_subs = {
x: r * sin(theta) * cos(phi),
y: r * sin(theta) * sin(phi),
z: r * cos(theta)
}
g_sph = g.subs(forward_subs).simplify()
$\left[\begin{matrix}1 & 0 & 0\\0 & \frac{1}{r^{2}} & 0\\0 & 0 & \frac{1}{r^{2} \sin^{2}{\left(\theta \right)}}\end{matrix}\right]$
This gives me the expected metric tensor as a matrix but there's a lot of extra effort:
forward_subs
since the CoordinateSymbol
"x" does not exist until I create Car3D
, and I have to provide the forward and backward transformations to create Car3D
.I can convert the metric tensor Matrix
to a tensor expression as follows:
g_sph_tens = 0
for i, e_i in enumerate(Car3D.base_oneforms()):
for j, e_j in enumerate(Car3D.base_oneforms()):
g_sph_tens += g_sph[i, j] * TP(e_i, e_j)
However, computing the Christoffel symbol for g_sph_tens
using metric_to_Christoffel_2nd
gives all zeros, which is incorrect.
Can somebody help me out with this? It seems like I'm missing some fundamentals.
Hi, I want to learn SymPy to do some calculations. I am practicing with the electromagnetic Lagrangian term $F_{\mu\nu}F^{\mu\nu}$. I have defined $F$ in the way that is explained in the documentation but cannot figure out how to go from there to $E^2 + B^2$. How am I supposed to do this?
Up to now I have this:
Ex, Ey, Ez, Bx, By, Bz = sp.symbols('E_x E_y E_z B_x B_y B_z')
c = sp.symbols('c', positive=True)
F = spT.TensorHead('F', [Lorentz, Lorentz], spT.TensorSymmetry.fully_symmetric(-2))
repl = {Lorentz: sp.diag(1, -1, -1, -1)}
repl.update({F(-mu,-nu): [
[0, Ex/c, Ey/c, Ez/c],
[-Ex/c, 0, -Bz, By],
[-Ey/c, Bz, 0, -Bx],
[-Ez/c, -By, Bx, 0]]}
)
and I have tried (F(mu,nu)*F(-mu,-nu)).replace_with_arrays(repl, [mu,nu])
and also F(mu,nu).replace_with_arrays(repl, [mu,nu])*F(-mu,-nu).replace_with_arrays(repl, [mu,nu])
.
sympy.preview(expr, output="pdf", viewer="file", filename="DFN_equations.pdf", dvioptions=['-D', '1100'], euler=False)
I wanted to expand on my question from earlier since it may be too specific.
I am working through some equations from the theory of elasticity.
The differential equation I am trying to work with is the equation
for compatibility of small displacements in a strained isotropic elastic body.
This is written as
$(1 - 2\eta) v^i|^j_j + v^j|^i_j = 0$
where
$f^i|_j$
represents the covariant derivative of f with respect to j, and
$f^i|^j$
represents the covariant derivative of f with respect to j
with the j index raised (this convention is from Green and Zerna -
I don't know if this is a standard notation).
We can use solutions of the above differential equation to compute
the stress tensor as follows:
$\frac{\tau^{ij}}{\mu} = g^{js} v^i|_s + g^{ir} v^j|_r + \frac{2 \eta}{1 - 2 \eta} g^{ij} v^r|_r$
One of the solutions to the displacement equation is as follows:
$v_i = F|_i$
where F is a harmonic function.
What's a good way to
\limits
is.
\limits
puts them directly above or below the integral sign. That does seem better for non-inline mode, especially if the limits can be long.
I am working in a simple QFT calculation and would like to do it using Sympy to learn (and also check my result). I have found the Quantum Mechanics module but cannot see how to start using it for my purpose. I have defined these quantities:
import sympy.physics.quantum as Q
vacuum = Q.OrthogonalKet(0)
annihilation_op = Q.Operator('a')
creation_op = Q.Dagger(annihilation_op)
and now I want to tell Sympy that $a |n\rangle = \sqrt{n} | n-1 \rangle$ if $n>0$ else $0$ and $a^\dagger | n \rangle = \sqrt{n+1} | n+1 \rangle$. How would I do this? Also, how do I impose the commutation relations between $a$ and $a^\dagger$?