print("Hello World!")
3**123
# Boolean XOR
True^False
# Division versus integer division in Python 3.
print(1/3)
print(1//3)
for w in ["bim", "bum", "bam"]:
word = ""
if w == "bim":
word = word + " BIM BIM"
elif w == "bum":
word += " BUM BUM"
else:
word += " BAM BAM !!!"
word = w + " " + word
print(word)
def factors(n):
"""
An inefficient but simple method for finding the factors of
an integer number n, omitting n itself.
"""
return [i for i in range(1, n//2+1) if n % i == 0]
factors(144)
# Help on that function
?factors
def perfects(n):
""" Computing perfect numbers less than n, using a previously defined function. """
return [i for i in range(1, n) if sum(factors(i)) == i]
perfects(2000)
# Start by importing all default objects from sympy
from sympy import *
#Set up best available printing interface
init_printing()
I, pi, E, EulerGamma, oo, sqrt(-2), log(-2), S(1)/3, Rational(1, 3)
beta = symbols('beta')
z = symbols('z_0:10')
beta, z
sin(z[3]) + beta**5
By default symbols are complex valued. If we do not want them to behave as complex numbers, we should specify another behaviour using the keywords:
commutative, complex, imaginary, real, integer, odd, even, prime, composite, zero, nonzero, rational, algebraic, transcendental, irrational, finite, infinite, negative, nonnegative, positive, nonpositive, hermitian, antihermitian
x = symbols('x')
n = symbols('n', integer = True)
sin(pi*x) + cos(pi*n)
Series, limits, derivatives, integrals.
x, y, z, t = symbols('x, y, z, t')
a, b, c = symbols('a, b, c')
m, n = symbols('m, n', integer = True)
[
summation(m**3, (m, 1, n)),
summation(1/n**17, (n, 1, oo)),
limit(sin(x)/x, x, 0),
limit(sin(exp(-1/x)+x)/x, x, 0, '+'),
diff(sin(x)**2, x),
diff(exp(x*y), x, 2, y),
integrate(cos(x), x),
integrate(exp(-x**2), (x, -oo, +oo))
]
""" A table of integrals. """
from IPython.display import Math, display
iData=[(sin(x)**2, (x, 0, pi)),
(t**(z-1)*exp(-t), (t, 0, oo)),
(exp(-t*x)*log(t), (t, 0, oo)),
((1-t)**(S(1)/2)*t**(S(1)/3), (t, 0, 1))]
for d in iData:
i = Integral(*d)
result = latex(i) + " = " + latex(simplify(i.doit()))
display(Math(result))
f1 = log(gamma(x))
f1s = series(f1, x, 0, 8)
result = (latex(f1) + r"\\=" + latex(f1s) + r"\\="
+ latex(f1s.rewrite(zeta)))
display(Math(result))
print('\n', result)
f2 = sin(x**Rational(1, 3) + x**Rational(1, 5))
f2s = series(f2, x, 0, 1.5)
display(Math(latex(f2) + "=" + latex(f2s)))
N(pi**E, 1000)
Let us sum the following series: $$\sum_{k=1}^{\infty}\left(\frac{1}{k}-\log(1+\frac{1}{k})\right) = \gamma_E$$
k = symbols('k', integer = True)
s1 = N(Sum(1/k - log(1 + 1/k), (k, 1, oo)), 20)
s2 = N(EulerGamma, 50)
print(s1)
print(s2)
print("%.3g" %(s1-s2))
An inverse Mellin tranform:
\begin{align} \int_{-\infty}^\infty dt\; \Gamma\left(it+\tfrac{1}{2}\right) \pi^{-it} = 2\pi^\frac{3}{2} e^{-\pi} \end{align}i1 = N(Integral(gamma(I*t + S(1)/2)*pi**(-I*t), (t, -oo, +oo)), 20)
i2 = N(2*pi**(S(3)/2)*exp(-pi), 30)
print(i1)
print(i2)
print("%.3g" %abs(i1-i2))
%matplotlib inline
plot(*[legendre(i, x) for i in range(1, 5)], (x, -1, 1))
from sympy.plotting.plot import *
r, phi = symbols('r, phi')
a = 20.1
b = 0.05
c = 10
d = 0.35
plot3d_parametric_surface(r*sin(a*r)*(1 + d*r*sin(phi)),
r*cos(a*r)*(1 + d*r*cos(phi)),
b*r**c,
(r, 0, 1), (phi, 0, 2*pi),
nb_of_points_u = 100,
nb_of_points_v = 100,
title = "Snake")
Let's do some object-oriented programming.
In 2d CFT, the Kac table is the set of dimension $$ \Delta_{(r,s)} = \frac14\left( \frac{q}{p}(r^2-1) +\frac{p}{q}(s^2-1) -2rs\right) \quad \text{with} \quad \left\{\begin{array}{l} 1\leq r\leq p-1 \\ 1\leq s\leq q-1 \end{array}\right. $$ for $p, q$ positive integers.
class Kac:
def __init__(self, indices = (3, 4)):
self.indices = indices
(p, q) = indices
self.table = [[self.Dimension((r, q-s)) for r in range(1, p)]
for s in range(1, q)]
def Dimension(self, pair):
(p, q) = self.indices
(r, s) = pair
return (S(q)/p*(r**2-1) + S(p)/q*(s**2-1) + 2*(1-r*s))/4
def display(self, source = False):
if source:
print(latex(Matrix(self.table)))
else:
return Matrix(self.table)
def ground_state(self):
return min(sum(self.table, []))
myKac = Kac((5, 8))
myKac.Dimension((3, 2))
myKac.display()
myKac.display(source = True)
myKac.ground_state()
from IPython.display import Math,Latex,display
from sympy import *
from sympy.diffgeom import *
# Let's define our own abbreviations for some useful functions
TP = TensorProduct
d = Differential
LieD = LieDerivative
# A manifold, a patch, and coordinates, with our own names
m = Manifold('Schwarzschild_4', 4)
p = Patch('External', m)
spherical = CoordSystem('Spherical', p, ['t', 'r', 'theta', 'phi'])
# Coordinate, derivatives (vector fields), one-forms
xs = t, r, theta, phi = spherical.coord_functions()
es = e_t, e_r, e_theta, e_phi = spherical.base_vectors()
fs = dt, dr, dtheta, dphi = spherical.base_oneforms()
xs, es, fs
Let's check that $dx^j(\partial_{x^i}) = \delta_i^j$
pprint(Matrix([[f(e) for f in fs] for e in es]))
We introduce a metric that depends on two unknown functions: $$ ds^2 = -A_0(r) dt^2 + A_1(r) dr^2 + r^2(d\theta^2 + \sin^2\theta d\phi^2) $$
(A0, A1) = symbols('A_0, A_1', cls = Function)
metric = (-A0(r)*TP(dt, dt) + A1(r)*TP(dr, dr)
+ r**2*(TP(dtheta, dtheta) + sin(theta)**2*TP(dphi, dphi)))
metric
""" We compute the Ricci tensor, check that its off-diagonal terms
identically vanish, and write the vanishing of the diagonal terms
as equations for A_0, A_1. """
ricci = metric_to_Ricci_components(metric)
print("Off-diagonal terms of the Ricci tensor:",
[[ricci[i, j] for i in range(4) if not(i == j)] for j in range(4)])
# Technicality: r is a ScalarField, for solving ODE's it is better
# to replace it by a symbol
R = symbols('R', Positive = True)
ricciII = ricci.replace(r, R).doit()
print("\nDiagonal terms of the Ricci tensor:")
for i in range(4):
display(Math("(" + str(i) + ")\quad "
+ latex(Eq(ricciII[i, i], 0))))
""" We find that the first two equations have a simple combination
that can be easily solved. """
eq5 = expand( R * A1(R) * (ricciII[0, 0]*A1(R) + ricciII[1, 1]*A0(R)) )
display(Eq(eq5, 0))
sol5 = dsolve(eq5, A1(R), ics = {A1(1):1/A0(1)})
# Specifying initial conditions
display(sol5)
""" We plug this solution into the original four equations. """
ricciIII = ricciII.replace(sol5.lhs, sol5.rhs).doit()
for i in range(4):
display(Math("(" + str(i) + "')\quad "
+ latex(Eq(ricciIII[i,i], 0))))
""" We solve (2'), and check that all equations are satisfied. """
M = Symbol('M')
sol = dsolve(ricciIII[2,2], A0(R), ics = {A0(M):0})
ricciIV = ricciIII.replace(sol.lhs, sol.rhs).doit()
display(sol)
for i in range(4):
display(Math("(" + str(i) + "')\quad "
+ latex(Eq(simplify(ricciIV[i, i]), 0))))
""" Computing the metric. """
schwMetric = ( metric.replace(r, R).replace(sol5.lhs, sol5.rhs).doit()
.replace(sol.lhs, sol.rhs).doit() )
schwMetric
SymPy: symbolic computing in Python Cite this PeerJ article if you use Sympy.