Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Browse files

Coding style: examples

  • Loading branch information...
commit b94ae1a286d5af7ed81a77a5e23effdc06511e49 1 parent 477fc2a
@jrioux jrioux authored
View
23 examples/advanced/autowrap_integrators.py
@@ -36,7 +36,7 @@
from sympy.utilities.lambdify import implemented_function
from sympy.utilities.autowrap import autowrap, ufuncify
from sympy import Idx, IndexedBase, Lambda, pprint, Symbol, oo, Integral,\
- Function
+ Function
from sympy.physics.sho import R_nl
from sympy.physics.hydrogen import R_nl as hydro_nl
@@ -78,19 +78,16 @@ def main():
expr = expr.evalf(15)
print "The h.o. wave function with l = %i and n = %i is" % (
- orbital_momentum_l, n)
+ orbital_momentum_l, n)
pprint(expr)
# implement, compile and wrap it as a ufunc
basis_ho[n] = ufuncify(x, expr)
-
# now let's see if we can express a hydrogen radial wave in terms of
# the ho basis. Here's the solution we will approximate:
-
H_ufunc = ufuncify(x, hydro_nl(hydrogen_n, orbital_momentum_l, 1, x))
-
# The transformation to a different basis can be written like this,
#
# psi(r) = sum_i c(i) phi_i(r)
@@ -110,7 +107,6 @@ def main():
# the low-level code. (In fact, summations are very easy to create, and as
# we will see it is often necessary to take extra steps in order to avoid
# them.)
-
# we need one integration ufunc for each wave function in the h.o. basis
binary_integrator = {}
for n in range(basis_dimension):
@@ -131,7 +127,6 @@ def main():
psi_ho = implemented_function('psi_ho',
Lambda(x, R_nl(n, orbital_momentum_l, omega2, x)))
-
# We represent the hydrogen function by an array which will be an input
# argument to the binary routine. This will let the integrators find
# h.o. basis coefficients for any wave function we throw at them.
@@ -150,7 +145,7 @@ def main():
expr = A[i]**2*psi_ho(A[i])*psi[i]*step
- if n==0:
+ if n == 0:
print "Setting up binary integrators for the integral:"
pprint(Integral(x**2*psi_ho(x)*Function('psi')(x), (x, 0, oo)))
@@ -165,28 +160,26 @@ def main():
binary_integrator[n] = autowrap(expr, args=[A.label, psi.label, step, m])
# Lets see how it converges with the grid dimension
- print "Checking convergence of integrator for n = %i" %n
+ print "Checking convergence of integrator for n = %i" % n
for g in range(3, 8):
grid, step = np.linspace(0, rmax, 2**g, retstep=True)
print "grid dimension %5i, integral = %e" % (2**g,
binary_integrator[n](grid, H_ufunc(grid), step))
-
print "A binary integrator has been set up for each basis state"
print "We will now use them to reconstruct a hydrogen solution."
# Note: We didn't need to specify grid or use gridsize before now
grid, stepsize = np.linspace(0, rmax, gridsize, retstep=True)
- print "Calculating coefficients with gridsize = %i and stepsize %f" %(
- len(grid), stepsize)
+ print "Calculating coefficients with gridsize = %i and stepsize %f" % (
+ len(grid), stepsize)
coeffs = {}
for n in range(basis_dimension):
coeffs[n] = binary_integrator[n](grid, H_ufunc(grid), stepsize)
print "c(%i) = %e" % (n, coeffs[n])
-
print "Constructing the approximate hydrogen wave"
hydro_approx = 0
all_steps = {}
@@ -194,8 +187,7 @@ def main():
hydro_approx += basis_ho[n](grid)*coeffs[n]
all_steps[n] = hydro_approx.copy()
if pylab:
- line = pylab.plot(grid, all_steps[n], ':', label='max n = %i'%n)
-
+ line = pylab.plot(grid, all_steps[n], ':', label='max n = %i' % n)
# check error numerically
diff = np.max(np.abs(hydro_approx - H_ufunc(grid)))
@@ -205,7 +197,6 @@ def main():
else:
print "Ah, that's a pretty good approximation!"
-
# Check visually
if pylab:
print "Here's a plot showing the contribution for each n"
View
12 examples/advanced/curvilinear_coordinates.py
@@ -10,7 +10,7 @@
"""
from sympy import var, sin, cos, pprint, Matrix, eye, trigsimp, Eq, \
- Function, simplify, sinh, cosh, expand
+ Function, simplify, sinh, cosh, expand
def laplace(f, g_inv, g_det, X):
"""
@@ -85,15 +85,15 @@ def main():
Matrix([rho*sin(theta)*cos(phi), rho*sin(theta)*sin(phi),
rho*cos(theta)]),
[rho, theta, phi],
- recursive=True
- )
+ recursive=True
+ )
transform("rotating disk",
- Matrix([t, x*cos(w*t)-y*sin(w*t), x*sin(w*t)+y*cos(w*t), z]),
+ 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]),
+ Matrix([sigma*tau, (tau**2 - sigma**2)/2]),
[sigma, tau])
# too complex:
@@ -106,7 +106,7 @@ def main():
transform("elliptic",
Matrix([a*cosh(mu)*cos(nu), a*sinh(mu)*sin(nu)]),
[mu, nu]
- )
+ )
if __name__ == "__main__":
main()
View
23 examples/advanced/dense_coding_example.py
@@ -10,7 +10,8 @@
from sympy.physics.quantum.grover import superposition_basis
def main():
- psi = superposition_basis(2); psi
+ psi = superposition_basis(2)
+ psi
# Dense coding demo:
@@ -25,26 +26,30 @@ def main():
# To Send Bob the message |0>|0>
print "To Send Bob the message |00>."
- circuit = H(1)*CNOT(1,0)
- result = qapply(circuit*psi); result
+ circuit = H(1)*CNOT(1, 0)
+ result = qapply(circuit*psi)
+ result
pprint(result)
# To send Bob the message |0>|1>
print "To Send Bob the message |01>."
- circuit = H(1)*CNOT(1,0)*X(1)
- result = qapply(circuit*psi); result
+ circuit = H(1)*CNOT(1, 0)*X(1)
+ result = qapply(circuit*psi)
+ result
pprint(result)
# To send Bob the message |1>|0>
print "To Send Bob the message |10>."
- circuit = H(1)*CNOT(1,0)*Z(1)
- result = qapply(circuit*psi); result
+ circuit = H(1)*CNOT(1, 0)*Z(1)
+ result = qapply(circuit*psi)
+ result
pprint(result)
# To send Bob the message |1>|1>
print "To Send Bob the message |11>."
- circuit = H(1)*CNOT(1,0)*Z(1)*X(1)
- result = qapply(circuit*psi); result
+ circuit = H(1)*CNOT(1, 0)*Z(1)*X(1)
+ result = qapply(circuit*psi)
+ result
pprint(result)
if __name__ == "__main__":
View
316 examples/advanced/fem.py
@@ -16,191 +16,187 @@
"""
from sympy import symbols, Symbol, factorial, Rational, zeros, div, eye, \
- integrate, diff, pprint, reduced
+ integrate, diff, pprint, reduced
x, y, z = symbols('x,y,z')
class ReferenceSimplex:
- def __init__(self, nsd):
- self.nsd = nsd
- if nsd <= 3:
- coords = symbols('x,y,z')[:nsd]
- else:
- coords = [Symbol("x_%d" % d) for d in range(nsd)]
- self.coords = coords
-
- def integrate(self,f):
- coords = self.coords
- nsd = self.nsd
-
- limit = 1
- for p in coords:
- limit -= p
-
- intf = f
- for d in range(0,nsd):
- p = coords[d]
- limit += p
- intf = integrate(intf, (p, 0, limit))
- return intf
+ def __init__(self, nsd):
+ self.nsd = nsd
+ if nsd <= 3:
+ coords = symbols('x,y,z')[:nsd]
+ else:
+ coords = [Symbol("x_%d" % d) for d in range(nsd)]
+ self.coords = coords
+
+ def integrate(self, f):
+ coords = self.coords
+ nsd = self.nsd
+
+ limit = 1
+ for p in coords:
+ limit -= p
+
+ intf = f
+ for d in range(0, nsd):
+ p = coords[d]
+ limit += p
+ intf = integrate(intf, (p, 0, limit))
+ return intf
def bernstein_space(order, nsd):
- if nsd > 3:
- raise RuntimeError("Bernstein only implemented in 1D, 2D, and 3D")
- sum = 0
- basis = []
- coeff = []
-
- if nsd == 1:
- b1, b2 = x, 1-x
- for o1 in range(0,order+1):
- for o2 in range(0,order+1):
- if o1 + o2 == order:
- aij = Symbol("a_%d_%d" % (o1,o2))
- sum += aij*binomial(order,o1)*pow(b1, o1)*pow(b2, o2)
- basis.append(binomial(order,o1)*pow(b1, o1)*pow(b2, o2))
- coeff.append(aij)
-
-
- if nsd == 2:
- b1, b2, b3 = x, y, 1-x-y
- for o1 in range(0,order+1):
- for o2 in range(0,order+1):
- for o3 in range(0,order+1):
- if o1 + o2 + o3 == order:
- aij = Symbol("a_%d_%d_%d" % (o1,o2,o3))
- fac = factorial(order) / (factorial(o1)*factorial(o2)*factorial(o3))
- sum += aij*fac*pow(b1, o1)*pow(b2, o2)*pow(b3, o3)
- basis.append(fac*pow(b1, o1)*pow(b2, o2)*pow(b3, o3))
- coeff.append(aij)
-
- if nsd == 3:
- b1, b2, b3, b4 = x, y, z, 1-x-y-z
- for o1 in range(0,order+1):
- for o2 in range(0,order+1):
- for o3 in range(0,order+1):
- for o4 in range(0,order+1):
- if o1 + o2 + o3 + o4 == order:
- aij = Symbol("a_%d_%d_%d_%d" % (o1,o2,o3,o4))
- fac = factorial(order)/ (factorial(o1)*factorial(o2)*factorial(o3)*factorial(o4))
- sum += aij*fac*pow(b1, o1)*pow(b2, o2)*pow(b3, o3)*pow(b4, o4)
- basis.append(fac*pow(b1, o1)*pow(b2, o2)*pow(b3, o3)*pow(b4, o4))
- coeff.append(aij)
-
-
- return sum, coeff, basis
+ if nsd > 3:
+ raise RuntimeError("Bernstein only implemented in 1D, 2D, and 3D")
+ sum = 0
+ basis = []
+ coeff = []
+
+ if nsd == 1:
+ b1, b2 = x, 1 - x
+ for o1 in range(0, order + 1):
+ for o2 in range(0, order + 1):
+ if o1 + o2 == order:
+ aij = Symbol("a_%d_%d" % (o1, o2))
+ sum += aij*binomial(order, o1)*pow(b1, o1)*pow(b2, o2)
+ basis.append(binomial(order, o1)*pow(b1, o1)*pow(b2, o2))
+ coeff.append(aij)
+
+ if nsd == 2:
+ b1, b2, b3 = x, y, 1 - x - y
+ for o1 in range(0, order + 1):
+ for o2 in range(0, order + 1):
+ for o3 in range(0, order + 1):
+ if o1 + o2 + o3 == order:
+ aij = Symbol("a_%d_%d_%d" % (o1, o2, o3))
+ fac = factorial(order) / (factorial(o1)*factorial(o2)*factorial(o3))
+ sum += aij*fac*pow(b1, o1)*pow(b2, o2)*pow(b3, o3)
+ basis.append(fac*pow(b1, o1)*pow(b2, o2)*pow(b3, o3))
+ coeff.append(aij)
+
+ if nsd == 3:
+ b1, b2, b3, b4 = x, y, z, 1 - x - y - z
+ for o1 in range(0, order + 1):
+ for o2 in range(0, order + 1):
+ for o3 in range(0, order + 1):
+ for o4 in range(0, order + 1):
+ if o1 + o2 + o3 + o4 == order:
+ aij = Symbol("a_%d_%d_%d_%d" % (o1, o2, o3, o4))
+ fac = factorial(order)/ (factorial(o1)*factorial(o2)*factorial(o3)*factorial(o4))
+ sum += aij*fac*pow(b1, o1)*pow(b2, o2)*pow(b3, o3)*pow(b4, o4)
+ basis.append(fac*pow(b1, o1)*pow(b2, o2)*pow(b3, o3)*pow(b4, o4))
+ coeff.append(aij)
+
+ return sum, coeff, basis
def create_point_set(order, nsd):
- h = Rational(1,order)
- set = []
-
- if nsd == 1:
- for i in range(0, order+1):
- x = i*h
- if x <= 1:
- set.append((x,y))
-
- if nsd == 2:
- for i in range(0, order+1):
- x = i*h
- for j in range(0, order+1):
- y = j*h
- if x + y <= 1:
- set.append((x,y))
-
- if nsd == 3:
- for i in range(0, order+1):
- x = i*h
- for j in range(0, order+1):
- y = j*h
- for k in range(0, order+1):
- z = j*h
- if x + y + z <= 1:
- set.append((x,y,z))
-
- return set
-
+ h = Rational(1, order)
+ set = []
+
+ if nsd == 1:
+ for i in range(0, order + 1):
+ x = i*h
+ if x <= 1:
+ set.append((x, y))
+
+ if nsd == 2:
+ for i in range(0, order + 1):
+ x = i*h
+ for j in range(0, order + 1):
+ y = j*h
+ if x + y <= 1:
+ set.append((x, y))
+
+ if nsd == 3:
+ for i in range(0, order + 1):
+ x = i*h
+ for j in range(0, order + 1):
+ y = j*h
+ for k in range(0, order + 1):
+ z = j*h
+ if x + y + z <= 1:
+ set.append((x, y, z))
+
+ return set
def create_matrix(equations, coeffs):
- A = zeros(len(equations))
- i = 0; j = 0
- for j in range(0, len(coeffs)):
- c = coeffs[j]
- for i in range(0, len(equations)):
- e = equations[i]
- d, _ = reduced(e, [c])
- A[i,j] = d[0]
- return A
-
+ A = zeros(len(equations))
+ i = 0
+ j = 0
+ for j in range(0, len(coeffs)):
+ c = coeffs[j]
+ for i in range(0, len(equations)):
+ e = equations[i]
+ d, _ = reduced(e, [c])
+ A[i, j] = d[0]
+ return A
class Lagrange:
- def __init__(self,nsd, order):
- self.nsd = nsd
- self.order = order
- self.compute_basis()
+ def __init__(self, nsd, order):
+ self.nsd = nsd
+ self.order = order
+ self.compute_basis()
- def nbf(self):
- return len(self.N)
+ def nbf(self):
+ return len(self.N)
- def compute_basis(self):
- order = self.order
- nsd = self.nsd
- N = []
- pol, coeffs, basis = bernstein_space(order, nsd)
- points = create_point_set(order, nsd)
+ def compute_basis(self):
+ order = self.order
+ nsd = self.nsd
+ N = []
+ pol, coeffs, basis = bernstein_space(order, nsd)
+ points = create_point_set(order, nsd)
- equations = []
- for p in points:
- ex = pol.subs(x, p[0])
- if nsd > 1:
- ex = ex.subs(y, p[1])
- if nsd > 2:
- ex = ex.subs(z, p[2])
- equations.append(ex )
+ equations = []
+ for p in points:
+ ex = pol.subs(x, p[0])
+ if nsd > 1:
+ ex = ex.subs(y, p[1])
+ if nsd > 2:
+ ex = ex.subs(z, p[2])
+ equations.append(ex )
- A = create_matrix(equations, coeffs)
- Ainv = A.inv()
+ A = create_matrix(equations, coeffs)
+ Ainv = A.inv()
- b = eye(len(equations))
+ b = eye(len(equations))
- xx = Ainv*b
+ xx = Ainv*b
- for i in range(0,len(equations)):
- Ni = pol
- for j in range(0,len(coeffs)):
- Ni = Ni.subs(coeffs[j], xx[j,i])
- N.append(Ni)
+ for i in range(0, len(equations)):
+ Ni = pol
+ for j in range(0, len(coeffs)):
+ Ni = Ni.subs(coeffs[j], xx[j, i])
+ N.append(Ni)
- self.N = N
+ self.N = N
def main():
- t = ReferenceSimplex(2)
- fe = Lagrange(2,2)
-
- u = 0
- #compute u = sum_i u_i N_i
- us = []
- for i in range(0, fe.nbf()):
- ui = Symbol("u_%d" % i)
- us.append(ui)
- u += ui*fe.N[i]
-
-
- J = zeros(fe.nbf())
- for i in range(0, fe.nbf()):
- Fi = u*fe.N[i]
- print Fi
- for j in range(0, fe.nbf()):
- uj = us[j]
- integrands = diff(Fi, uj)
- print integrands
- J[j,i] = t.integrate(integrands)
-
- pprint(J)
+ t = ReferenceSimplex(2)
+ fe = Lagrange(2, 2)
+
+ u = 0
+ #compute u = sum_i u_i N_i
+ us = []
+ for i in range(0, fe.nbf()):
+ ui = Symbol("u_%d" % i)
+ us.append(ui)
+ u += ui*fe.N[i]
+
+ J = zeros(fe.nbf())
+ for i in range(0, fe.nbf()):
+ Fi = u*fe.N[i]
+ print Fi
+ for j in range(0, fe.nbf()):
+ uj = us[j]
+ integrands = diff(Fi, uj)
+ print integrands
+ J[j, i] = t.integrate(integrands)
+
+ pprint(J)
if __name__ == "__main__":
- main()
+ main()
View
12 examples/advanced/gibbs_phenomenon.py
@@ -14,7 +14,7 @@
"""
from sympy import var, sqrt, integrate, conjugate, seterr, Abs, pprint, I, pi,\
- sin, cos, sign, Plot, lambdify, Integral, S
+ sin, cos, sign, Plot, lambdify, Integral, S
#seterr(True)
@@ -52,7 +52,7 @@ def l2_projection(f, basis, lim):
"""
r = 0
for b in basis:
- r += l2_inner_product(f, b, lim) * b
+ r += l2_inner_product(f, b, lim) * b
return r
def l2_gram_schmidt(list, lim):
@@ -104,14 +104,14 @@ def msolve(f, x):
f = lambdify(x, f)
x0 = -0.001
dx = 0.001
- while f(x0-dx) * f(x0) > 0:
- x0 = x0-dx
- x_max = x0-dx
+ while f(x0 - dx) * f(x0) > 0:
+ x0 = x0 - dx
+ x_max = x0 - dx
x_min = x0
assert f(x_max) > 0
assert f(x_min) < 0
for n in range(100):
- x0 = (x_max+x_min)/2
+ x0 = (x_max + x_min)/2
if f(x0) > 0:
x_max = x0
else:
View
16 examples/advanced/pidigits.py
@@ -17,20 +17,20 @@ def display_fraction(digits, skip=0, colwidth=10, columns=5):
"""Pretty printer for first n digits of a fraction"""
perline = colwidth * columns
printed = 0
- for linecount in range((len(digits)-skip) // (colwidth * columns)):
- line = digits[skip+linecount*perline:skip+(linecount+1)*perline]
+ for linecount in range((len(digits) - skip) // (colwidth * columns)):
+ line = digits[skip + linecount*perline:skip + (linecount + 1)*perline]
for i in range(columns):
- print line[i*colwidth : (i+1)*colwidth],
- print ":", (linecount+1)*perline
- if (linecount+1) % 10 == 0:
+ print line[i*colwidth: (i + 1)*colwidth],
+ print ":", (linecount + 1)*perline
+ if (linecount + 1) % 10 == 0:
print
printed += colwidth*columns
- rem = (len(digits)-skip) % (colwidth * columns)
+ rem = (len(digits) - skip) % (colwidth * columns)
if rem:
buf = digits[-rem:]
s = ""
for i in range(columns):
- s += buf[:colwidth].ljust(colwidth+1, " ")
+ s += buf[:colwidth].ljust(colwidth + 1, " ")
buf = buf[colwidth:]
print s + ":", printed + colwidth*columns
@@ -43,7 +43,7 @@ def calculateit(func, base, n, tofile):
else:
skip = len(intpart)
print "Step 1 of 2: calculating binary value..."
- prec = int(n*math.log(base,2))+10
+ prec = int(n*math.log(base, 2)) + 10
t = clock()
a = func(prec)
step1_time = clock() - t
View
72 examples/advanced/pyglet_plotting.py
@@ -15,7 +15,7 @@
from time import sleep, clock
def main():
- x,y,z = symbols('x,y,z')
+ x, y, z = symbols('x,y,z')
# toggle axes visibility with F5, colors with F6
axes_options = 'visible=false; colored=true; label_ticks=true; label_axes=true; overlay=true; stride=0.5'
@@ -24,19 +24,20 @@ def main():
p = PygletPlot(width=600, height=500, ortho=False, invert_mouse_zoom=False, axes=axes_options, antialiasing=True)
examples = []
+
def example_wrapper(f):
examples.append(f)
return f
@example_wrapper
def mirrored_saddles():
- p[5] = x**2-y**2, [20], [20]
- p[6] = y**2-x**2, [20], [20]
+ p[5] = x**2 - y**2, [20], [20]
+ p[6] = y**2 - x**2, [20], [20]
@example_wrapper
def mirrored_saddles_saveimage():
- p[5] = x**2-y**2, [20], [20]
- p[6] = y**2-x**2, [20], [20]
+ p[5] = x**2 - y**2, [20], [20]
+ p[6] = y**2 - x**2, [20], [20]
p.wait_for_calculations()
# although the calculation is complete,
# we still need to wait for it to be
@@ -46,19 +47,19 @@ def mirrored_saddles_saveimage():
@example_wrapper
def mirrored_ellipsoids():
- p[2] = x**2+y**2, [40], [40], 'color=zfade'
- p[3] = -x**2-y**2, [40], [40], 'color=zfade'
+ p[2] = x**2 + y**2, [40], [40], 'color=zfade'
+ p[3] = -x**2 - y**2, [40], [40], 'color=zfade'
@example_wrapper
def saddle_colored_by_derivative():
- f = x**2-y**2
+ f = x**2 - y**2
p[1] = f, 'style=solid'
p[1].color = abs(f.diff(x)), abs(f.diff(x) + f.diff(y)), abs(f.diff(y))
@example_wrapper
def ding_dong_surface():
- f = sqrt(1.0-y)*y
- p[1] = f, [x,0,2*pi,40], [y,-1,4,100], 'mode=cylindrical; style=solid; color=zfade4'
+ f = sqrt(1.0 - y)*y
+ p[1] = f, [x, 0, 2*pi, 40], [y, -1, 4, 100], 'mode=cylindrical; style=solid; color=zfade4'
@example_wrapper
def polar_circle():
@@ -67,7 +68,7 @@ def polar_circle():
@example_wrapper
def polar_flower():
p[8] = 1.5*sin(4*x), [160], 'mode=polar'
- p[8].color = z, x, y, (0.5,0.5,0.5), (0.8,0.8,0.8), (x,y,None,z) # z is used for t
+ p[8].color = z, x, y, (0.5, 0.5, 0.5), (0.8, 0.8, 0.8), (x, y, None, z) # z is used for t
@example_wrapper
def simple_cylinder():
@@ -76,39 +77,39 @@ def simple_cylinder():
@example_wrapper
def cylindrical_hyperbola():
## (note that polar is an alias for cylindrical)
- p[10] = 1/y, 'mode=polar', [x], [y,-2,2,20]
+ p[10] = 1/y, 'mode=polar', [x], [y, -2, 2, 20]
@example_wrapper
def extruded_hyperbolas():
- p[11] = 1/x, [x,-10,10,100], [1], 'style=solid'
- p[12] = -1/x, [x,-10,10,100], [1], 'style=solid'
+ p[11] = 1/x, [x, -10, 10, 100], [1], 'style=solid'
+ p[12] = -1/x, [x, -10, 10, 100], [1], 'style=solid'
@example_wrapper
def torus():
- a,b = 1, 0.5 # radius, thickness
- p[13] = (a+b*cos(x))*cos(y), (a+b*cos(x))*sin(y), b*sin(x), [x,0,pi*2,40], [y,0,pi*2,40]
+ a, b = 1, 0.5 # radius, thickness
+ p[13] = (a + b*cos(x))*cos(y), (a + b*cos(x))*sin(y), b*sin(x), [x, 0, pi*2, 40], [y, 0, pi*2, 40]
@example_wrapper
def warped_torus():
- a,b = 2, 1 # radius, thickness
- p[13] = (a+b*cos(x))*cos(y), (a+b*cos(x))*sin(y), b*sin(x)+0.5*sin(4*y), [x,0,pi*2,40], [y,0,pi*2,40]
+ a, b = 2, 1 # radius, thickness
+ p[13] = (a + b*cos(x))*cos(y), (a + b*cos(x))*sin(y), b*sin(x) + 0.5*sin(4*y), [x, 0, pi*2, 40], [y, 0, pi*2, 40]
@example_wrapper
def parametric_spiral():
- p[14] = cos(y), sin(y), y/10.0, [y,-4*pi,4*pi,100]
- p[14].color = x,(0.1,0.9),y,(0.1,0.9),z,(0.1,0.9)
+ p[14] = cos(y), sin(y), y/10.0, [y, -4*pi, 4*pi, 100]
+ p[14].color = x, (0.1, 0.9), y, (0.1, 0.9), z, (0.1, 0.9)
@example_wrapper
def multistep_gradient():
p[1] = 1, 'mode=spherical', 'style=both'
#p[1] = exp(-x**2-y**2+(x*y)/4), [-1.7,1.7,100], [-1.7,1.7,100], 'style=solid'
#p[1] = 5*x*y*exp(-x**2-y**2), [-2,2,100], [-2,2,100]
- gradient = [ 0.0, (0.3, 0.3, 1.0),
+ gradient = [ 0.0, (0.3, 0.3, 1.0),
0.30, (0.3, 1.0, 0.3),
- 0.55, (0.95,1.0, 0.2),
- 0.65, (1.0,0.95, 0.2),
+ 0.55, (0.95, 1.0, 0.2),
+ 0.65, (1.0, 0.95, 0.2),
0.85, (1.0, 0.7, 0.2),
- 1.0, (1.0, 0.3, 0.2) ]
+ 1.0, (1.0, 0.3, 0.2) ]
p[1].color = z, [None, None, z], gradient
#p[1].color = 'zfade'
#p[1].color = 'zfade3'
@@ -116,14 +117,14 @@ def multistep_gradient():
@example_wrapper
def lambda_vs_sympy_evaluation():
start = clock()
- p[4] = x**2+y**2, [100], [100], 'style=solid'
+ p[4] = x**2 + y**2, [100], [100], 'style=solid'
p.wait_for_calculations()
- print "lambda-based calculation took %s seconds." % (clock()-start)
+ print "lambda-based calculation took %s seconds." % (clock() - start)
start = clock()
- p[4] = x**2+y**2, [100], [100], 'style=solid; use_sympy_eval'
+ p[4] = x**2 + y**2, [100], [100], 'style=solid; use_sympy_eval'
p.wait_for_calculations()
- print "sympy substitution-based calculation took %s seconds." % (clock()-start)
+ print "sympy substitution-based calculation took %s seconds." % (clock() - start)
@example_wrapper
def gradient_vectors():
@@ -139,8 +140,8 @@ def draw_gradient_vectors(f, iu, iv):
representing the gradient of f.
"""
dx, dy, dz = f.diff(x), f.diff(y), 0
- FF = lambdify( [x,y], [x,y,f] )
- FG = lambdify( [x,y], [dx,dy,dz] )
+ FF = lambdify( [x, y], [x, y, f] )
+ FG = lambdify( [x, y], [dx, dy, dz] )
iu.v_steps /= 5
iv.v_steps /= 5
Gvl = list(list([FF(u, v), FG(u, v)]
@@ -171,7 +172,7 @@ def draw():
glEnd()
return draw
- p[i] = f, [-0.5,0.5,25], [-0.5,0.5,25], 'style=solid'
+ p[i] = f, [-0.5, 0.5, 25], [-0.5, 0.5, 25], 'style=solid'
iu = PlotInterval(p[i].intervals[0])
iv = PlotInterval(p[i].intervals[1])
p[i].postdraw.append(draw_gradient_vectors(f, iu, iv))
@@ -180,9 +181,9 @@ def draw():
gradient_vectors_inner(-x**2 - y**2, 2)
def help_str():
- s = ("\nPlot p has been created. Useful commands: \n"
- " help(p), p[1] = x**2, print p, p.clear() \n\n"
- "Available examples (see source in plotting.py):\n\n")
+ s = ("\nPlot p has been created. Useful commands: \n"
+ " help(p), p[1] = x**2, print p, p.clear() \n\n"
+ "Available examples (see source in plotting.py):\n\n")
for i in xrange(len(examples)):
s += "(%i) %s\n" % (i, examples[i].__name__)
s += "\n"
@@ -197,7 +198,8 @@ def example(i):
elif i >= 0 and i < len(examples):
p.clear()
examples[i]()
- else: print "Not a valid example.\n"
+ else:
+ print "Not a valid example.\n"
print p
example(0) # 0 - 15 are defined above
View
101 examples/advanced/qft.py
@@ -15,56 +15,56 @@
"""
-from sympy import Basic,exp,Symbol,sin,Rational,I,Mul, Matrix, \
+from sympy import Basic, exp, Symbol, sin, Rational, I, Mul, Matrix, \
ones, sqrt, pprint, simplify, Eq, sympify
from sympy.physics import msigma, mgamma
#gamma^mu
-gamma0=mgamma(0)
-gamma1=mgamma(1)
-gamma2=mgamma(2)
-gamma3=mgamma(3)
-gamma5=mgamma(5)
+gamma0 = mgamma(0)
+gamma1 = mgamma(1)
+gamma2 = mgamma(2)
+gamma3 = mgamma(3)
+gamma5 = mgamma(5)
#sigma_i
-sigma1=msigma(1)
-sigma2=msigma(2)
-sigma3=msigma(3)
+sigma1 = msigma(1)
+sigma2 = msigma(2)
+sigma3 = msigma(3)
E = Symbol("E", real=True)
m = Symbol("m", real=True)
-def u(p,r):
+def u(p, r):
""" p = (p1, p2, p3); r = 0,1 """
- assert r in [1,2]
- p1,p2,p3 = p
+ assert r in [1, 2]
+ p1, p2, p3 = p
if r == 1:
- ksi = Matrix([[1],[0]])
+ ksi = Matrix([[1], [0]])
else:
- ksi = Matrix([[0],[1]])
- a = (sigma1*p1 + sigma2*p2 + sigma3*p3) / (E+m) * ksi
- if a ==0:
+ ksi = Matrix([[0], [1]])
+ a = (sigma1*p1 + sigma2*p2 + sigma3*p3) / (E + m) * ksi
+ if a == 0:
a = zeros(2, 1)
- return sqrt(E+m) * Matrix([[ksi[0,0]], [ksi[1,0]], [a[0,0]], [a[1,0]]])
+ return sqrt(E + m) * Matrix([[ksi[0, 0]], [ksi[1, 0]], [a[0, 0]], [a[1, 0]]])
-def v(p,r):
+def v(p, r):
""" p = (p1, p2, p3); r = 0,1 """
- assert r in [1,2]
- p1,p2,p3 = p
+ assert r in [1, 2]
+ p1, p2, p3 = p
if r == 1:
- ksi = Matrix([[1],[0]])
+ ksi = Matrix([[1], [0]])
else:
- ksi = -Matrix([[0],[1]])
- a = (sigma1*p1 + sigma2*p2 + sigma3*p3) / (E+m) * ksi
- if a ==0:
+ ksi = -Matrix([[0], [1]])
+ a = (sigma1*p1 + sigma2*p2 + sigma3*p3) / (E + m) * ksi
+ if a == 0:
a = zeros(2, 1)
- return sqrt(E+m) * Matrix([[a[0,0]], [a[1,0]], [ksi[0,0]], [ksi[1,0]]])
+ return sqrt(E + m) * Matrix([[a[0, 0]], [a[1, 0]], [ksi[0, 0]], [ksi[1, 0]]])
def pslash(p):
- p1,p2,p3 = p
- p0 = sqrt(m**2+p1**2+p2**2+p3**2)
- return gamma0*p0-gamma1*p1-gamma2*p2-gamma3*p3
+ p1, p2, p3 = p
+ p0 = sqrt(m**2 + p1**2 + p2**2 + p3**2)
+ return gamma0*p0 - gamma1*p1 - gamma2*p2 - gamma3*p3
def Tr(M):
return M.trace()
@@ -73,56 +73,55 @@ def xprint(lhs, rhs):
pprint( Eq(sympify(lhs), rhs ) )
def main():
- a=Symbol("a", real=True)
- b=Symbol("b", real=True)
- c=Symbol("c", real=True)
+ a = Symbol("a", real=True)
+ b = Symbol("b", real=True)
+ c = Symbol("c", real=True)
- p = (a,b,c)
+ p = (a, b, c)
assert u(p, 1).D * u(p, 2) == Matrix(1, 1, [0])
assert u(p, 2).D * u(p, 1) == Matrix(1, 1, [0])
- p1,p2,p3 =[Symbol(x, real=True) for x in ["p1","p2","p3"]]
- pp1,pp2,pp3 =[Symbol(x, real=True) for x in ["pp1","pp2","pp3"]]
- k1,k2,k3 =[Symbol(x, real=True) for x in ["k1","k2","k3"]]
- kp1,kp2,kp3 =[Symbol(x, real=True) for x in ["kp1","kp2","kp3"]]
+ p1, p2, p3 = [Symbol(x, real=True) for x in ["p1", "p2", "p3"]]
+ pp1, pp2, pp3 = [Symbol(x, real=True) for x in ["pp1", "pp2", "pp3"]]
+ k1, k2, k3 = [Symbol(x, real=True) for x in ["k1", "k2", "k3"]]
+ kp1, kp2, kp3 = [Symbol(x, real=True) for x in ["kp1", "kp2", "kp3"]]
- p = (p1,p2,p3)
- pp = (pp1,pp2,pp3)
+ p = (p1, p2, p3)
+ pp = (pp1, pp2, pp3)
- k = (k1,k2,k3)
- kp = (kp1,kp2,kp3)
+ k = (k1, k2, k3)
+ kp = (kp1, kp2, kp3)
mu = Symbol("mu")
- e = (pslash(p)+m*ones(4))*(pslash(k)-m*ones(4))
- f = pslash(p)+m*ones(4)
- g = pslash(p)-m*ones(4)
-
+ e = (pslash(p) + m*ones(4))*(pslash(k) - m*ones(4))
+ f = pslash(p) + m*ones(4)
+ g = pslash(p) - m*ones(4)
#pprint(e)
xprint( 'Tr(f*g)', Tr(f*g) )
#print Tr(pslash(p) * pslash(k)).expand()
- M0 = [ ( v(pp, 1).D * mgamma(mu) * u(p, 1) ) * ( u(k, 1).D * mgamma(mu,True) * \
+ M0 = [ ( v(pp, 1).D * mgamma(mu) * u(p, 1) ) * ( u(k, 1).D * mgamma(mu, True) *
v(kp, 1) ) for mu in range(4)]
- M = M0[0]+M0[1]+M0[2]+M0[3]
+ M = M0[0] + M0[1] + M0[2] + M0[3]
M = M[0]
assert isinstance(M, Basic)
#print M
#print simplify(M)
- d=Symbol("d", real=True) #d=E+m
+ d = Symbol("d", real=True) # d=E+m
xprint('M', M)
print "-"*40
- M = ((M.subs(E,d-m)).expand() * d**2 ).expand()
- xprint('M2', 1/(E+m)**2 * M)
+ M = ((M.subs(E, d - m)).expand() * d**2 ).expand()
+ xprint('M2', 1/(E + m)**2 * M)
print "-"*40
- x,y= M.as_real_imag()
+ x, y = M.as_real_imag()
xprint('Re(M)', x)
xprint('Im(M)', y)
- e = x**2+y**2
+ e = x**2 + y**2
xprint('abs(M)**2', e)
print "-"*40
xprint('Expand(abs(M)**2)', e.expand())
View
157 examples/advanced/relativity.py
@@ -16,79 +16,79 @@
from sympy import (exp, Symbol, sin, Rational, Derivative, dsolve, Function,
Matrix, Eq, pprint, Pow, classify_ode, solve)
-def grad(f,X):
- a=[]
+def grad(f, X):
+ a = []
for x in X:
a.append(f.diff(x))
return a
-def d(m,x):
- return grad(m[0,0],x)
+def d(m, x):
+ return grad(m[0, 0], x)
class MT(object):
- def __init__(self,m):
- self.gdd=m
- self.guu=m.inv()
+ def __init__(self, m):
+ self.gdd = m
+ self.guu = m.inv()
def __str__(self):
return "g_dd =\n" + str(self.gdd)
- def dd(self,i,j):
- return self.gdd[i,j]
+ def dd(self, i, j):
+ return self.gdd[i, j]
- def uu(self,i,j):
- return self.guu[i,j]
+ def uu(self, i, j):
+ return self.guu[i, j]
class G(object):
- def __init__(self,g,x):
+ def __init__(self, g, x):
self.g = g
self.x = x
- def udd(self,i,k,l):
- g=self.g
- x=self.x
- r=0
- for m in [0,1,2,3]:
- r+=g.uu(i,m)/2 * (g.dd(m,k).diff(x[l])+g.dd(m,l).diff(x[k]) \
- - g.dd(k,l).diff(x[m]))
+ def udd(self, i, k, l):
+ g = self.g
+ x = self.x
+ r = 0
+ for m in [0, 1, 2, 3]:
+ r += g.uu(i, m)/2 * (g.dd(m, k).diff(x[l]) + g.dd(m, l).diff(x[k])
+ - g.dd(k, l).diff(x[m]))
return r
class Riemann(object):
- def __init__(self,G,x):
+ def __init__(self, G, x):
self.G = G
self.x = x
- def uddd(self,rho,sigma,mu,nu):
- G=self.G
- x=self.x
- r=G.udd(rho,nu,sigma).diff(x[mu])-G.udd(rho,mu,sigma).diff(x[nu])
- for lam in [0,1,2,3]:
- r+=G.udd(rho,mu,lam)*G.udd(lam,nu,sigma) \
- -G.udd(rho,nu,lam)*G.udd(lam,mu,sigma)
+ def uddd(self, rho, sigma, mu, nu):
+ G = self.G
+ x = self.x
+ r = G.udd(rho, nu, sigma).diff(x[mu]) - G.udd(rho, mu, sigma).diff(x[nu])
+ for lam in [0, 1, 2, 3]:
+ r += G.udd(rho, mu, lam)*G.udd(lam, nu, sigma) \
+ - G.udd(rho, nu, lam)*G.udd(lam, mu, sigma)
return r
class Ricci(object):
- def __init__(self,R,x):
+ def __init__(self, R, x):
self.R = R
self.x = x
self.g = R.G.g
- def dd(self,mu,nu):
- R=self.R
- x=self.x
- r=0
- for lam in [0,1,2,3]:
- r+=R.uddd(lam,mu,lam,nu)
+ def dd(self, mu, nu):
+ R = self.R
+ x = self.x
+ r = 0
+ for lam in [0, 1, 2, 3]:
+ r += R.uddd(lam, mu, lam, nu)
return r
- def ud(self,mu,nu):
- r=0
- for lam in [0,1,2,3]:
- r+=self.g.uu(mu,lam)*self.dd(lam,nu)
+ def ud(self, mu, nu):
+ r = 0
+ for lam in [0, 1, 2, 3]:
+ r += self.g.uu(mu, lam)*self.dd(lam, nu)
return r.expand()
def curvature(Rmn):
- return Rmn.ud(0,0)+Rmn.ud(1,1)+Rmn.ud(2,2)+Rmn.ud(3,3)
+ return Rmn.ud(0, 0) + Rmn.ud(1, 1) + Rmn.ud(2, 2) + Rmn.ud(3, 3)
#class nu(Function):
# def getname(self):
@@ -102,18 +102,18 @@ def curvature(Rmn):
nu = Function("nu")
lam = Function("lambda")
-t=Symbol("t")
-r=Symbol("r")
-theta=Symbol(r"theta")
-phi=Symbol(r"phi")
+t = Symbol("t")
+r = Symbol("r")
+theta = Symbol(r"theta")
+phi = Symbol(r"phi")
#general, spherically symmetric metric
-gdd=Matrix((
- (-exp(nu(r)),0,0,0),
+gdd = Matrix((
+ (-exp(nu(r)), 0, 0, 0),
(0, exp(lam(r)), 0, 0),
(0, 0, r**2, 0),
(0, 0, 0, r**2*sin(theta)**2)
- ))
+))
#spherical - flat
#gdd=Matrix((
# (-1, 0, 0, 0),
@@ -135,86 +135,85 @@ def curvature(Rmn):
# (0, 0, r**2*sin(theta)**2, 0),
# (0, 0, 0, r**2)
# ))
-g=MT(gdd)
-X=(t,r,theta,phi)
-Gamma=G(g,X)
-Rmn=Ricci(Riemann(Gamma,X),X)
+g = MT(gdd)
+X = (t, r, theta, phi)
+Gamma = G(g, X)
+Rmn = Ricci(Riemann(Gamma, X), X)
-def pprint_Gamma_udd(i,k,l):
- pprint(Eq(Symbol('Gamma^%i_%i%i' % (i,k,l)), Gamma.udd(i,k,l)))
+def pprint_Gamma_udd(i, k, l):
+ pprint(Eq(Symbol('Gamma^%i_%i%i' % (i, k, l)), Gamma.udd(i, k, l)))
-def pprint_Rmn_dd(i,j):
- pprint(Eq(Symbol('R_%i%i' % (i,j)), Rmn.dd(i,j)))
+def pprint_Rmn_dd(i, j):
+ pprint(Eq(Symbol('R_%i%i' % (i, j)), Rmn.dd(i, j)))
# from Differential Equations example
def eq1():
r = Symbol("r")
- e = Rmn.dd(0,0)
+ e = Rmn.dd(0, 0)
e = e.subs(nu(r), -lam(r))
pprint(dsolve(e, lam(r)))
def eq2():
r = Symbol("r")
- e = Rmn.dd(1,1)
+ e = Rmn.dd(1, 1)
C = Symbol("CC")
e = e.subs(nu(r), -lam(r))
pprint(dsolve(e, lam(r)))
def eq3():
r = Symbol("r")
- e = Rmn.dd(2,2)
+ e = Rmn.dd(2, 2)
e = e.subs(nu(r), -lam(r))
pprint(dsolve(e, lam(r)))
def eq4():
r = Symbol("r")
- e = Rmn.dd(3,3)
+ e = Rmn.dd(3, 3)
e = e.subs(nu(r), -lam(r))
pprint(dsolve(e, lam(r)))
pprint(dsolve(e, lam(r), 'best'))
-
def main():
print "Initial metric:"
pprint(gdd)
print "-"*40
print "Christoffel symbols:"
- pprint_Gamma_udd(0,1,0)
- pprint_Gamma_udd(0,0,1)
+ pprint_Gamma_udd(0, 1, 0)
+ pprint_Gamma_udd(0, 0, 1)
print
- pprint_Gamma_udd(1,0,0)
- pprint_Gamma_udd(1,1,1)
- pprint_Gamma_udd(1,2,2)
- pprint_Gamma_udd(1,3,3)
+ pprint_Gamma_udd(1, 0, 0)
+ pprint_Gamma_udd(1, 1, 1)
+ pprint_Gamma_udd(1, 2, 2)
+ pprint_Gamma_udd(1, 3, 3)
print
- pprint_Gamma_udd(2,2,1)
- pprint_Gamma_udd(2,1,2)
- pprint_Gamma_udd(2,3,3)
+ pprint_Gamma_udd(2, 2, 1)
+ pprint_Gamma_udd(2, 1, 2)
+ pprint_Gamma_udd(2, 3, 3)
print
- pprint_Gamma_udd(3,2,3)
- pprint_Gamma_udd(3,3,2)
- pprint_Gamma_udd(3,1,3)
- pprint_Gamma_udd(3,3,1)
+ pprint_Gamma_udd(3, 2, 3)
+ pprint_Gamma_udd(3, 3, 2)
+ pprint_Gamma_udd(3, 1, 3)
+ pprint_Gamma_udd(3, 3, 1)
print"-"*40
print"Ricci tensor:"
- pprint_Rmn_dd(0,0)
- e = Rmn.dd(1,1)
- pprint_Rmn_dd(1,1)
- pprint_Rmn_dd(2,2)
- pprint_Rmn_dd(3,3)
+ pprint_Rmn_dd(0, 0)
+ e = Rmn.dd(1, 1)
+ pprint_Rmn_dd(1, 1)
+ pprint_Rmn_dd(2, 2)
+ pprint_Rmn_dd(3, 3)
#print
#print "scalar curvature:"
#print curvature(Rmn)
print "-"*40
print "Solve Einstein's equations:"
e = e.subs(nu(r), -lam(r)).doit()
- l = dsolve(e, lam(r))
+ l = dsolve(e, lam(r))
pprint(l)
lamsol = solve(l, lam(r))[0]
- metric = gdd.subs(lam(r), lamsol).subs(nu(r),-lamsol)#.combine()
+ metric = gdd.subs(lam(r), lamsol).subs(nu(r), -lamsol) # .combine()
print "metric:"
pprint(metric)
View
7 examples/all.py
@@ -71,7 +71,7 @@
"advanced.pidigits",
"advanced.qft",
"advanced.relativity",
- ]
+]
WINDOWED_EXAMPLES = [
"beginner.plotting_nice_plot",
@@ -81,7 +81,7 @@
"advanced.autowrap_integrators",
"advanced.autowrap_ufuncify",
"advanced.pyglet_plotting",
- ]
+]
EXAMPLE_DIR = os.path.dirname(__file__)
@@ -177,7 +177,8 @@ def run_example(example, reporter=None):
class DummyFile(object):
- def write(self, x): pass
+ def write(self, x):
+ pass
def suppress_output(fn):
View
16 examples/beginner/limits_examples.py
@@ -18,21 +18,21 @@ def main():
a = Symbol("a")
h = Symbol("h")
- show( limit(sqrt(x**2 - 5*x + 6) - x, x, oo) , -Rational(5)/2 )
+ show( limit(sqrt(x**2 - 5*x + 6) - x, x, oo), -Rational(5)/2 )
- show( limit(x*(sqrt(x**2 + 1) - x), x, oo) , Rational(1)/2 )
+ show( limit(x*(sqrt(x**2 + 1) - x), x, oo), Rational(1)/2 )
- show( limit(x - sqrt3(x**3 - 1), x, oo) , Rational(0) )
+ show( limit(x - sqrt3(x**3 - 1), x, oo), Rational(0) )
- show( limit(log(1 + exp(x))/x, x, -oo) , Rational(0) )
+ show( limit(log(1 + exp(x))/x, x, -oo), Rational(0) )
- show( limit(log(1 + exp(x))/x, x, oo) , Rational(1) )
+ show( limit(log(1 + exp(x))/x, x, oo), Rational(1) )
- show( limit(sin(3*x)/x, x, 0) , Rational(3) )
+ show( limit(sin(3*x)/x, x, 0), Rational(3) )
- show( limit(sin(5*x)/sin(2*x), x, 0) , Rational(5)/2 )
+ show( limit(sin(5*x)/sin(2*x), x, 0), Rational(5)/2 )
- show( limit(((x - 1)/(x + 1))**x, x, oo) , exp(-2))
+ show( limit(((x - 1)/(x + 1))**x, x, oo), exp(-2))
if __name__ == "__main__":
main()
View
22 examples/beginner/plot_examples.py
@@ -3,7 +3,7 @@
from sympy import Symbol, exp, sin, cos
from sympy.plotting import (plot, plot_parametric,
- plot3d_parametric_surface,plot3d_parametric_line,
+ plot3d_parametric_surface, plot3d_parametric_line,
plot3d)
lx = range(5)
@@ -16,28 +16,28 @@
expr = x**2 - 1
b = plot(expr, (x, 2, 4), show=False) # cartesian plot
-e = plot(exp(-x),(x, 0, 4), show=False) # cartesian plot (and coloring, see below)
+e = plot(exp(-x), (x, 0, 4), show=False) # cartesian plot (and coloring, see below)
f = plot3d_parametric_line(sin(x), cos(x), x, (x, 0, 10), show=False) # 3d parametric line plot
g = plot3d(sin(x)*cos(y), (x, -5, 5), (y, -10, 10), show=False) # 3d surface cartesian plot
-h = plot3d_parametric_surface(cos(u)*v, sin(u)*v, u, (u, 0, 10),(v, -2, 2), show=False) # 3d parametric surface plot
+h = plot3d_parametric_surface(cos(u)*v, sin(u)*v, u, (u, 0, 10), (v, -2, 2), show=False) # 3d parametric surface plot
# Some aesthetics
-e[0].line_color = lambda x : x / 4
-f[0].line_color = lambda x, y, z : z / 10
-g[0].surface_color = lambda x, y : sin(x)
+e[0].line_color = lambda x: x / 4
+f[0].line_color = lambda x, y, z: z / 10
+g[0].surface_color = lambda x, y: sin(x)
# Some more stuff on aesthetics - coloring wrt coordinates or parameters
param_line_2d = plot_parametric((x*cos(x), x*sin(x), (x, 0, 15)), (1.1*x*cos(x), 1.1*x*sin(x), (x, 0, 15)), show=False)
-param_line_2d[0].line_color = lambda u : sin(u) # parametric
-param_line_2d[1].line_color = lambda u, v : u**2 + v**2 # coordinates
+param_line_2d[0].line_color = lambda u: sin(u) # parametric
+param_line_2d[1].line_color = lambda u, v: u**2 + v**2 # coordinates
param_line_2d.title = 'The inner one is colored by parameter and the outher one by coordinates'
param_line_3d = plot3d_parametric_line((x*cos(x), x*sin(x), x, (x, 0, 15)),
(1.5*x*cos(x), 1.5*x*sin(x), x, (x, 0, 15)),
(2*x*cos(x), 2*x*sin(x), x, (x, 0, 15)), show=False)
-param_line_3d[0].line_color = lambda u : u #parametric
-param_line_3d[1].line_color = lambda u, v : u*v # first and second coordinates
-param_line_3d[2].line_color = lambda u, v, w : u*v*w # all coordinates
+param_line_3d[0].line_color = lambda u: u # parametric
+param_line_3d[1].line_color = lambda u, v: u*v # first and second coordinates
+param_line_3d[2].line_color = lambda u, v, w: u*v*w # all coordinates
if __name__ == '__main__':
View
64 examples/intermediate/coupled_cluster.py
@@ -14,29 +14,29 @@
symbols, expand, pprint, Rational, latex, Dummy
)
-pretty_dummies_dict={
- 'above':'cdefgh',
- 'below':'klmno',
- 'general':'pqrstu'
- }
+pretty_dummies_dict = {
+ 'above': 'cdefgh',
+ 'below': 'klmno',
+ 'general': 'pqrstu'
+}
def get_CC_operators():
"""
Returns a tuple (T1,T2) of unique operators.
"""
- i = symbols('i',below_fermi=True,cls=Dummy)
- a = symbols('a',above_fermi=True,cls=Dummy)
- t_ai = AntiSymmetricTensor('t',(a,),(i,))
+ i = symbols('i', below_fermi=True, cls=Dummy)
+ a = symbols('a', above_fermi=True, cls=Dummy)
+ t_ai = AntiSymmetricTensor('t', (a,), (i,))
ai = NO(Fd(a)*F(i))
- i,j = symbols('i,j',below_fermi=True,cls=Dummy)
- a,b = symbols('a,b',above_fermi=True,cls=Dummy)
- t_abij = AntiSymmetricTensor('t',(a,b),(i,j))
+ i, j = symbols('i,j', below_fermi=True, cls=Dummy)
+ a, b = symbols('a,b', above_fermi=True, cls=Dummy)
+ t_abij = AntiSymmetricTensor('t', (a, b), (i, j))
abji = NO(Fd(a)*Fd(b)*F(j)*F(i))
T1 = t_ai*ai
T2 = Rational(1, 4)*t_abij*abji
- return (T1,T2)
+ return (T1, T2)
def main():
print
@@ -47,10 +47,10 @@ def main():
print
# setup hamiltonian
- p,q,r,s = symbols('p,q,r,s',cls=Dummy)
- f = AntiSymmetricTensor('f',(p,),(q,))
+ p, q, r, s = symbols('p,q,r,s', cls=Dummy)
+ f = AntiSymmetricTensor('f', (p,), (q,))
pr = NO((Fd(p)*F(q)))
- v = AntiSymmetricTensor('v',(p,q),(r,s))
+ v = AntiSymmetricTensor('v', (p, q), (r, s))
pqsr = NO(Fd(p)*Fd(q)*F(s)*F(r))
H = f*pr + Rational(1, 4)*v*pqsr
@@ -59,36 +59,36 @@ def main():
print "Calculating 4 nested commutators"
C = Commutator
- T1,T2 = get_CC_operators()
- T = T1+ T2
+ T1, T2 = get_CC_operators()
+ T = T1 + T2
print "commutator 1..."
- comm1 = wicks(C(H,T))
+ comm1 = wicks(C(H, T))
comm1 = evaluate_deltas(comm1)
comm1 = substitute_dummies(comm1)
- T1,T2 = get_CC_operators()
- T = T1+ T2
+ T1, T2 = get_CC_operators()
+ T = T1 + T2
print "commutator 2..."
- comm2 = wicks(C(comm1,T))
+ comm2 = wicks(C(comm1, T))
comm2 = evaluate_deltas(comm2)
comm2 = substitute_dummies(comm2)
- T1,T2 = get_CC_operators()
- T = T1+ T2
+ T1, T2 = get_CC_operators()
+ T = T1 + T2
print "commutator 3..."
- comm3 = wicks(C(comm2,T))
+ comm3 = wicks(C(comm2, T))
comm3 = evaluate_deltas(comm3)
comm3 = substitute_dummies(comm3)
- T1,T2 = get_CC_operators()
- T = T1+ T2
+ T1, T2 = get_CC_operators()
+ T = T1 + T2
print "commutator 4..."
- comm4 = wicks(C(comm3,T))
+ comm4 = wicks(C(comm3, T))
comm4 = evaluate_deltas(comm4)
comm4 = substitute_dummies(comm4)
print "construct Hausdoff expansion..."
- eq = H + comm1+comm2/2 +comm3/6+comm4/24
+ eq = H + comm1 + comm2/2 + comm3/6 + comm4/24
eq = eq.expand()
eq = evaluate_deltas(eq)
eq = substitute_dummies(eq, new_indices=True,
@@ -97,8 +97,8 @@ def main():
print
print "extracting CC equations from full Hbar"
- i,j,k,l = symbols('i,j,k,l',below_fermi=True)
- a,b,c,d = symbols('a,b,c,d',above_fermi=True)
+ i, j, k, l = symbols('i,j,k,l', below_fermi=True)
+ a, b, c, d = symbols('a,b,c,d', above_fermi=True)
print
print "CC Energy:"
print latex(wicks(eq, simplify_dummies=True,
@@ -110,9 +110,9 @@ def main():
print latex(eqT1)
print
print "CC T2:"
- eqT2 = wicks(NO(Fd(i)*Fd(j)*F(b)*F(a))*eq,simplify_dummies=True, keep_only_fully_contracted=True, simplify_kronecker_deltas=True)
+ eqT2 = wicks(NO(Fd(i)*Fd(j)*F(b)*F(a))*eq, simplify_dummies=True, keep_only_fully_contracted=True, simplify_kronecker_deltas=True)
P = PermutationOperator
- eqT2 = simplify_index_permutations(eqT2,[P(a,b),P(i,j)])
+ eqT2 = simplify_index_permutations(eqT2, [P(a, b), P(i, j)])
print latex(eqT2)
if __name__ == "__main__":
View
4 examples/intermediate/infinite_1d_box.py
@@ -68,8 +68,8 @@ def energy_corrections(perturbation, n, a=10, mass=0.5):
Vnm(n, n, a).evalf(),
- (Vnm(n, n-1, a)**2/(E_n(n, a, mass) - E_n(n-1, a, mass))
- + Vnm(n, n+1, a)**2/(E_n(n, a, mass) - E_n(n+1, a, mass))).evalf())
+ (Vnm(n, n - 1, a)**2/(E_n(n, a, mass) - E_n(n - 1, a, mass))
+ + Vnm(n, n + 1, a)**2/(E_n(n, a, mass) - E_n(n + 1, a, mass))).evalf())
def main():
print
View
2  examples/intermediate/mplot2d.py
@@ -26,7 +26,7 @@ def mplot2d(f, var, show=True):
sys.exit("Matplotlib is required to use mplot2d.")
if not is_sequence(f):
- f = [f,]
+ f = [f, ]
for f_i in f:
x, y = sample(f_i, var)
View
6 examples/intermediate/mplot3d.py
@@ -24,7 +24,7 @@ def mplot3d(f, var1, var2, show=True):
p = import_module('pylab')
# Try newer version first
p3 = import_module('mpl_toolkits.mplot3d',
- __import__kwargs={'fromlist':['something']}) or import_module('matplotlib.axes3d')
+ __import__kwargs={'fromlist': ['something']}) or import_module('matplotlib.axes3d')
if not p or not p3:
sys.exit("Matplotlib is required to use mplot3d.")
@@ -34,7 +34,7 @@ def mplot3d(f, var1, var2, show=True):
ax = p3.Axes3D(fig)
#ax.plot_surface(x,y,z) #seems to be a bug in matplotlib
- ax.plot_wireframe(x,y,z)
+ ax.plot_wireframe(x, y, z)
ax.set_xlabel('X')
ax.set_ylabel('Y')
@@ -47,7 +47,7 @@ def main():
x = Symbol('x')
y = Symbol('y')
- mplot3d(x**2-y**2, (x, -10.0, 10.0, 20), (y, -10.0, 10.0, 20))
+ mplot3d(x**2 - y**2, (x, -10.0, 10.0, 20), (y, -10.0, 10.0, 20))
#mplot3d(x**2+y**2, (x, -10.0, 10.0, 20), (y, -10.0, 10.0, 20))
#mplot3d(sin(x)+sin(y), (x, -3.14, 3.14, 10), (y, -3.14, 3.14, 10))
View
48 examples/intermediate/partial_differential_eqs.py
@@ -14,45 +14,43 @@ def main():
R, Phi, Theta, u = map(Function, ['R', 'Phi', 'Theta', 'u'])
C1, C2 = symbols('C1,C2')
- pprint ("Separation of variables in Laplace equation in spherical coordinates")
- pprint ("Laplace equation in spherical coordinates:")
- eq = Eq(D(Xi(r, phi, theta), r, 2) + 2/r * D(Xi(r, phi, theta),r) + \
- 1/(r**2 * sin(phi)**2) * D(Xi(r, phi, theta), theta, 2) + \
- cos(phi)/(r**2 * sin(phi)) * D(Xi(r, phi, theta), phi) + \
+ pprint("Separation of variables in Laplace equation in spherical coordinates")
+ pprint("Laplace equation in spherical coordinates:")
+ eq = Eq(D(Xi(r, phi, theta), r, 2) + 2/r * D(Xi(r, phi, theta), r) +
+ 1/(r**2 * sin(phi)**2) * D(Xi(r, phi, theta), theta, 2) +
+ cos(phi)/(r**2 * sin(phi)) * D(Xi(r, phi, theta), phi) +
1/r**2 * D(Xi(r, phi, theta), phi, 2))
- pprint (eq)
+ pprint(eq)
- pprint ("We can either separate this equation in regards with variable r:")
+ pprint("We can either separate this equation in regards with variable r:")
res_r = pde_separate(eq, Xi(r, phi, theta), [R(r), u(phi, theta)])
- pprint (res_r)
+ pprint(res_r)
- pprint ("Or separate it in regards of theta:")
+ pprint("Or separate it in regards of theta:")
res_theta = pde_separate(eq, Xi(r, phi, theta), [Theta(theta), u(r, phi)])
- pprint (res_theta)
+ pprint(res_theta)
res_phi = pde_separate(eq, Xi(r, phi, theta), [Phi(phi), u(r, theta)])
- pprint ("But we cannot separate it in regards of variable phi: ")
- pprint ("Result: %s" % res_phi)
+ pprint("But we cannot separate it in regards of variable phi: ")
+ pprint("Result: %s" % res_phi)
- pprint ("\n\nSo let's make theta dependent part equal with -C1:")
+ pprint("\n\nSo let's make theta dependent part equal with -C1:")
eq_theta = Eq(res_theta[0], -C1)
- pprint (eq_theta)
- pprint ("\nThis also means that second part is also equal to -C1:")
+ pprint(eq_theta)
+ pprint("\nThis also means that second part is also equal to -C1:")
eq_left = Eq(res_theta[1], -C1)
- pprint (eq_left)
+ pprint(eq_left)
- pprint ("\nLets try to separate phi again :)")
+ pprint("\nLets try to separate phi again :)")
res_theta = pde_separate(eq_left, u(r, phi), [Phi(phi), R(r)])
- pprint ("\nThis time it is successful:")
- pprint (res_theta)
-
- pprint ("\n\nSo our final equations with separated variables are:")
- pprint (eq_theta)
- pprint (Eq(res_theta[0],C2))
- pprint (Eq(res_theta[1],C2))
-
+ pprint("\nThis time it is successful:")
+ pprint(res_theta)
+ pprint("\n\nSo our final equations with separated variables are:")
+ pprint(eq_theta)
+ pprint(Eq(res_theta[0], C2))
+ pprint(Eq(res_theta[1], C2))
if __name__ == "__main__":
View
6 examples/intermediate/sample.py
@@ -28,7 +28,7 @@ def sample2d(f, x_args):
x_l = float(x_max - x_min)
x_d = x_l/float(x_n)
- X = arange(float(x_min), float(x_max)+x_d, x_d)
+ X = arange(float(x_min), float(x_max) + x_d, x_d)
Y = empty(len(X))
for i in range(len(X)):
@@ -61,11 +61,11 @@ def sample3d(f, x_args, y_args):
x_l = float(x_max - x_min)
x_d = x_l/float(x_n)
- x_a = arange(float(x_min), float(x_max)+x_d, x_d)
+ x_a = arange(float(x_min), float(x_max) + x_d, x_d)
y_l = float(y_max - y_min)
y_d = y_l/float(y_n)
- y_a = arange(float(y_min), float(y_max)+y_d, y_d)
+ y_a = arange(float(y_min), float(y_max) + y_d, y_d)
def meshgrid(x, y):
"""
View
6 examples/intermediate/trees.py
@@ -15,14 +15,14 @@
from sympy import Symbol, Poly
def T(x):
- return x+x**2+2*x**3 + 4*x**4 + 9*x**5 + 20*x**6 + 48 * x**7 + \
- 115* x**8 + 286*x**9+719*x**10
+ return x + x**2 + 2*x**3 + 4*x**4 + 9*x**5 + 20*x**6 + 48 * x**7 + \
+ 115* x**8 + 286*x**9 + 719*x**10
def A(x):
return 1 + T(x) - T(x)**2/2 + T(x**2)/2
def main():
- x=Symbol("x")
+ x = Symbol("x")
s = Poly(A(x), x)
num = list(reversed(s.coeffs()))[:11]
View
22 examples/intermediate/vandermonde.py
@@ -63,7 +63,7 @@ def vandermonde(order, dim=1, syms='a b c d'):
for i in range(rank):
row_syms = [g.next() for g in generators]
all_syms.append(row_syms)
- for j,term in enumerate(terms):
+ for j, term in enumerate(terms):
v_entry = 1
for k in term:
v_entry *= row_syms[k]
@@ -82,8 +82,8 @@ def gen_poly(points, order, syms):
V, tmp_syms, terms = vandermonde(order, dim)
if num_pts < V.shape[0]:
raise ValueError(
- "Must provide %d points for order %d, dimension "\
- "%d polynomial, given %d points" % \
+ "Must provide %d points for order %d, dimension "
+ "%d polynomial, given %d points" %
(V.shape[0], order, dim, num_pts))
elif num_pts > V.shape[0]:
print "gen_poly given %d points but only requires %d, "\
@@ -101,7 +101,7 @@ def gen_poly(points, order, syms):
coeffs = V_inv.multiply(Matrix([points[i][-1] for i in xrange(num_pts)]))
f = 0
- for j,term in enumerate(terms):
+ for j, term in enumerate(terms):
t = 1
for k in term:
t *= syms[k]
@@ -128,21 +128,21 @@ def main():
\sum = %(sum)s
= %(sum_expand)s
""" % { "det": V.det(),
- "sum": det_sum,
- "sum_expand": det_sum.expand(),
+ "sum": det_sum,
+ "sum_expand": det_sum.expand(),
}
print '-'*79
print "Polynomial fitting with a Vandermonde Matrix:"
- x,y,z = symbols('x,y,z')
+ x, y, z = symbols('x,y,z')
- points = [(0,3), (1,2), (2,3)]
+ points = [(0, 3), (1, 2), (2, 3)]
print """
Quadratic function, represented by 3 points:
points = %(pts)s
f = %(f)s
""" % { "pts" : points,
- "f" : gen_poly(points, 2, [x]),
+ "f": gen_poly(points, 2, [x]),
}
points = [(0, 1, 1), (1, 0, 0), (1, 1, 0), (Rational(1, 2), 0, 0),
@@ -152,7 +152,7 @@ def main():
points = %(pts)s
f = %(f)s
""" % { "pts" : points,
- "f" : gen_poly(points, 2, [x, y]),
+ "f": gen_poly(points, 2, [x, y]),
}
points = [(0, 1, 1, 1), (1, 1, 0, 0), (1, 0, 1, 0), (1, 1, 1, 1)]
@@ -161,7 +161,7 @@ def main():
points = %(pts)s
f = %(f)s
""" % { "pts" : points,
- "f" : gen_poly(points, 1, [x, y, z]),
+ "f": gen_poly(points, 1, [x, y, z]),
}
Please sign in to comment.
Something went wrong with that request. Please try again.