def test_recursive():
from sympy import symbols, exp_polar, expand
a, b, c = symbols('a b c', positive=True)
r = exp(-(x - a)**2)*exp(-(x - b)**2)
e = integrate(r, (x, 0, oo), meijerg=True)
assert simplify(e.expand()) == (
sqrt(2)*sqrt(pi)*(
(erf(sqrt(2)*(a + b)/2) + 1)*exp(-a**2/2 + a*b - b**2/2))/4)
e = integrate(exp(-(x - a)**2)*exp(-(x - b)**2)*exp(c*x), (x, 0, oo), meijerg=True)
assert simplify(e) == (
sqrt(2)*sqrt(pi)*(erf(sqrt(2)*(2*a + 2*b + c)/4) + 1)*exp(-a**2 - b**2
+ (2*a + 2*b + c)**2/8)/4)
assert simplify(integrate(exp(-(x - a - b - c)**2), (x, 0, oo), meijerg=True)) == \
sqrt(pi)/2*(1 + erf(a + b + c))
assert simplify(integrate(exp(-(x + a + b + c)**2), (x, 0, oo), meijerg=True)) == \
sqrt(pi)/2*(1 - erf(a + b + c))
python类symbols()的实例源码
def runtest_ufuncify(language, backend):
has_module('numpy')
a, b, c = symbols('a b c')
fabc = ufuncify([a, b, c], a*b + c, language=language, backend=backend)
facb = ufuncify([a, c, b], a*b + c, language=language, backend=backend)
grid = numpy.linspace(-2, 2, 50)
for b in numpy.linspace(-5, 4, 3):
for c in numpy.linspace(-1, 1, 3):
expected = grid*b + c
assert numpy.sum(numpy.abs(expected - fabc(grid, b, c))) < 1e-13
assert numpy.sum(numpy.abs(expected - facb(grid, c, b))) < 1e-13
#
# tests of language-backend combinations
#
# f2py
def test_functions():
one_var = (acosh, ln, Heaviside, factorial, bernoulli, coth, tanh,
sign, arg, asin, DiracDelta, re, Abs, sinh, cos, cot, acos, acot,
gamma, bell, harmonic, LambertW, zeta, log, factorial, asinh,
acoth, cosh, dirichlet_eta, loggamma, erf, ceiling, im, fibonacci,
conjugate, tan, floor, atanh, sin, atan, lucas, exp)
two_var = (rf, ff, lowergamma, chebyshevu, chebyshevt, binomial,
atan2, polygamma, hermite, legendre, uppergamma)
x, y, z = symbols("x,y,z")
others = (chebyshevt_root, chebyshevu_root, Eijk(x, y, z),
Piecewise( (0, x < -1), (x**2, x <= 1), (x**3, True)),
assoc_legendre)
for cls in one_var:
check(cls)
c = cls(x)
check(c)
for cls in two_var:
check(cls)
c = cls(x, y)
check(c)
for cls in others:
check(cls)
#================== geometry ====================
def test_Routine_argument_order():
a, x, y, z = symbols('a x y z')
expr = (x + y)*z
raises(CodeGenArgumentListError, lambda: Routine("test", expr,
argument_sequence=[z, x]))
raises(CodeGenArgumentListError, lambda: Routine("test", Eq(a,
expr), argument_sequence=[z, x, y]))
r = Routine('test', Eq(a, expr), argument_sequence=[z, x, a, y])
assert [ arg.name for arg in r.arguments ] == [z, x, a, y]
assert [ type(arg) for arg in r.arguments ] == [
InputArgument, InputArgument, OutputArgument, InputArgument ]
r = Routine('test', Eq(z, expr), argument_sequence=[z, x, y])
assert [ type(arg) for arg in r.arguments ] == [
InOutArgument, InputArgument, InputArgument ]
from sympy.tensor import IndexedBase, Idx
A, B = map(IndexedBase, ['A', 'B'])
m = symbols('m', integer=True)
i = Idx('i', m)
r = Routine('test', Eq(A[i], B[i]), argument_sequence=[B, A, m])
assert [ arg.name for arg in r.arguments ] == [B.label, A.label, m]
expr = Integral(x*y*z, (x, 1, 2), (y, 1, 3))
r = Routine('test', Eq(a, expr), argument_sequence=[z, x, a, y])
assert [ arg.name for arg in r.arguments ] == [z, x, a, y]
def test_simple_c_codegen():
x, y, z = symbols('x,y,z')
expr = (x + y)*z
result = codegen(("test", expr), "C", "file", header=False, empty=False)
expected = [
("file.c",
"#include \"file.h\"\n"
"#include <math.h>\n"
"double test(double x, double y, double z) {\n"
" return z*(x + y);\n"
"}\n"),
("file.h",
"#ifndef PROJECT__FILE__H\n"
"#define PROJECT__FILE__H\n"
"double test(double x, double y, double z);\n"
"#endif\n")
]
assert result == expected
def test_dummy_loops_c():
from sympy.tensor import IndexedBase, Idx
# the following line could also be
# [Dummy(s, integer=True) for s in 'im']
# or [Dummy(integer=True) for s in 'im']
i, m = symbols('i m', integer=True, cls=Dummy)
x = IndexedBase('x')
y = IndexedBase('y')
i = Idx(i, m)
expected = (
'#include "file.h"\n'
'#include <math.h>\n'
'void test_dummies(int m_%(mno)i, double *x, double *y) {\n'
' for (int i_%(ino)i=0; i_%(ino)i<m_%(mno)i; i_%(ino)i++){\n'
' y[i_%(ino)i] = x[i_%(ino)i];\n'
' }\n'
'}\n'
) % {'ino': i.label.dummy_index, 'mno': m.dummy_index}
r = Routine('test_dummies', Eq(y[i], x[i]))
c = CCodeGen()
code = get_string(c.dump_c, [r])
assert code == expected
def test_output_arg_c():
from sympy import sin, cos, Equality
x, y, z = symbols("x,y,z")
r = Routine("foo", [Equality(y, sin(x)), cos(x)])
c = CCodeGen()
result = c.write([r], "test", header=False, empty=False)
assert result[0][0] == "test.c"
expected = (
'#include "test.h"\n'
'#include <math.h>\n'
'double foo(double x, double &y) {\n'
' y = sin(x);\n'
' return cos(x);\n'
'}\n'
)
assert result[0][1] == expected
def test_simple_f_code():
x, y, z = symbols('x,y,z')
expr = (x + y)*z
routine = Routine("test", expr)
code_gen = FCodeGen()
source = get_string(code_gen.dump_f95, [routine])
expected = (
"REAL*8 function test(x, y, z)\n"
"implicit none\n"
"REAL*8, intent(in) :: x\n"
"REAL*8, intent(in) :: y\n"
"REAL*8, intent(in) :: z\n"
"test = z*(x + y)\n"
"end function\n"
)
assert source == expected
def test_f_code_argument_order():
x, y, z = symbols('x,y,z')
expr = x + y
routine = Routine("test", expr, argument_sequence=[z, x, y])
code_gen = FCodeGen()
source = get_string(code_gen.dump_f95, [routine])
expected = (
"REAL*8 function test(z, x, y)\n"
"implicit none\n"
"REAL*8, intent(in) :: z\n"
"REAL*8, intent(in) :: x\n"
"REAL*8, intent(in) :: y\n"
"test = x + y\n"
"end function\n"
)
assert source == expected
def test_simple_f_header():
x, y, z = symbols('x,y,z')
expr = (x + y)*z
routine = Routine("test", expr)
code_gen = FCodeGen()
source = get_string(code_gen.dump_h, [routine])
expected = (
"interface\n"
"REAL*8 function test(x, y, z)\n"
"implicit none\n"
"REAL*8, intent(in) :: x\n"
"REAL*8, intent(in) :: y\n"
"REAL*8, intent(in) :: z\n"
"end function\n"
"end interface\n"
)
assert source == expected
def test_simple_f_codegen():
x, y, z = symbols('x,y,z')
expr = (x + y)*z
result = codegen(
("test", expr), "F95", "file", header=False, empty=False)
expected = [
("file.f90",
"REAL*8 function test(x, y, z)\n"
"implicit none\n"
"REAL*8, intent(in) :: x\n"
"REAL*8, intent(in) :: y\n"
"REAL*8, intent(in) :: z\n"
"test = z*(x + y)\n"
"end function\n"),
("file.h",
"interface\n"
"REAL*8 function test(x, y, z)\n"
"implicit none\n"
"REAL*8, intent(in) :: x\n"
"REAL*8, intent(in) :: y\n"
"REAL*8, intent(in) :: z\n"
"end function\n"
"end interface\n")
]
assert result == expected
def test_dummy_loops_f95():
from sympy.tensor import IndexedBase, Idx
# the following line could also be
# [Dummy(s, integer=True) for s in 'im']
# or [Dummy(integer=True) for s in 'im']
i, m = symbols('i m', integer=True, cls=Dummy)
x = IndexedBase('x')
y = IndexedBase('y')
i = Idx(i, m)
expected = (
'subroutine test_dummies(m_%(mcount)i, x, y)\n'
'implicit none\n'
'INTEGER*4, intent(in) :: m_%(mcount)i\n'
'REAL*8, intent(in), dimension(1:m_%(mcount)i) :: x\n'
'REAL*8, intent(out), dimension(1:m_%(mcount)i) :: y\n'
'INTEGER*4 :: i_%(icount)i\n'
'do i_%(icount)i = 1, m_%(mcount)i\n'
' y(i_%(icount)i) = x(i_%(icount)i)\n'
'end do\n'
'end subroutine\n'
) % {'icount': i.label.dummy_index, 'mcount': m.dummy_index}
r = Routine('test_dummies', Eq(y[i], x[i]))
c = FCodeGen()
code = get_string(c.dump_f95, [r])
assert code == expected
def test_output_arg_f():
from sympy import sin, cos, Equality
x, y, z = symbols("x,y,z")
r = Routine("foo", [Equality(y, sin(x)), cos(x)])
c = FCodeGen()
result = c.write([r], "test", header=False, empty=False)
assert result[0][0] == "test.f90"
assert result[0][1] == (
'REAL*8 function foo(x, y)\n'
'implicit none\n'
'REAL*8, intent(in) :: x\n'
'REAL*8, intent(out) :: y\n'
'y = sin(x)\n'
'foo = cos(x)\n'
'end function\n'
)
def Simple_manifold_with_scalar_function_derivative():
Print_Function()
coords = (x, y, z) = symbols('x y z')
basis = (e1, e2, e3, grad) = MV.setup('e_1 e_2 e_3', metric='[1,1,1]', coords=coords)
# Define surface
mfvar = (u, v) = symbols('u v')
X = u*e1 + v*e2 + (u**2 + v**2)*e3
print(X)
MF = Manifold(X, mfvar)
# Define field on the surface.
g = (v + 1)*log(u)
# Method 1: Using old Manifold routines.
VectorDerivative = (MF.rbasis[0]/MF.E_sq)*diff(g, u) + (MF.rbasis[1]/MF.E_sq)*diff(g, v)
print('Vector derivative =', VectorDerivative.subs({u: 1, v: 0}))
# Method 2: Using new Manifold routines.
dg = MF.Grad(g)
print('Vector derivative =', dg.subs({u: 1, v: 0}))
return
def derivatives_in_spherical_coordinates():
Print_Function()
X = (r, th, phi) = symbols('r theta phi')
curv = [[r*cos(phi)*sin(th), r*sin(phi)*sin(th), r*cos(th)], [1, r, r*sin(th)]]
(er, eth, ephi, grad) = MV.setup('e_r e_theta e_phi', metric='[1,1,1]', coords=X, curv=curv)
f = MV('f', 'scalar', fct=True)
A = MV('A', 'vector', fct=True)
B = MV('B', 'grade2', fct=True)
print('f =', f)
print('A =', A)
print('B =', B)
print('grad*f =', grad*f)
print('grad|A =', grad | A)
print('-I*(grad^A) =', -MV.I*(grad ^ A))
print('grad^B =', grad ^ B)
return
def Dirac_Equation_in_Geometric_Calculus():
Print_Function()
vars = symbols('t x y z')
(g0, g1, g2, g3, grad) = MV.setup('gamma*t|x|y|z', metric='[1,-1,-1,-1]', coords=vars)
I = MV.I
(m, e) = symbols('m e')
psi = MV('psi', 'spinor', fct=True)
A = MV('A', 'vector', fct=True)
sig_z = g3*g0
print('\\text{4-Vector Potential\\;\\;}\\bm{A} =', A)
print('\\text{8-component real spinor\\;\\;}\\bm{\\psi} =', psi)
dirac_eq = (grad*psi)*I*sig_z - e*A*psi - m*psi*g0
dirac_eq.simplify()
dirac_eq.Fmt(3, r'%\text{Dirac Equation\;\;}\nabla \bm{\psi} I \sigma_{z}-e\bm{A}\bm{\psi}-m\bm{\psi}\gamma_{t} = 0')
return
def derivatives_in_spherical_coordinates():
Print_Function()
X = (r, th, phi) = symbols('r theta phi')
curv = [[r*cos(phi)*sin(th), r*sin(phi)*sin(th), r*cos(th)], [1, r, r*sin(th)]]
(er, eth, ephi, grad) = MV.setup('e_r e_theta e_phi', metric='[1,1,1]', coords=X, curv=curv)
f = MV('f', 'scalar', fct=True)
A = MV('A', 'vector', fct=True)
B = MV('B', 'grade2', fct=True)
print('f =', f)
print('A =', A)
print('B =', B)
print('grad*f =', grad*f)
print('grad|A =', grad | A)
print('-I*(grad^A) =', (-MV.I*(grad ^ A)).simplify())
print('grad^B =', grad ^ B)
def MV_setup_options():
Print_Function()
(e1, e2, e3) = MV.setup('e_1 e_2 e_3', '[1,1,1]')
v = MV('v', 'vector')
print(v)
(e1, e2, e3) = MV.setup('e*1|2|3', '[1,1,1]')
v = MV('v', 'vector')
print(v)
(e1, e2, e3) = MV.setup('e*x|y|z', '[1,1,1]')
v = MV('v', 'vector')
print(v)
coords = symbols('x y z')
(e1, e2, e3, grad) = MV.setup('e', '[1,1,1]', coords=coords)
v = MV('v', 'vector')
print(v)
return
def main():
Format()
a = Matrix(2, 2, (1, 2, 3, 4))
b = Matrix(2, 1, (5, 6))
c = a*b
print(a, b, '=', c)
x, y = symbols('x, y')
d = Matrix(1, 2, (x**3, y**3))
e = Matrix(2, 2, (x**2, 2*x*y, 2*x*y, y**2))
f = d*e
print('%', d, e, '=', f)
xdvi()
return
def Simple_manifold_with_scalar_function_derivative():
Print_Function()
coords = (x, y, z) = symbols('x y z')
basis = (e1, e2, e3, grad) = MV.setup('e_1 e_2 e_3', metric='[1,1,1]', coords=coords)
# Define surface
mfvar = (u, v) = symbols('u v')
X = u*e1 + v*e2 + (u**2 + v**2)*e3
print('\\f{X}{u,v} =', X)
MF = Manifold(X, mfvar)
(eu, ev) = MF.Basis()
# Define field on the surface.
g = (v + 1)*log(u)
print('\\f{g}{u,v} =', g)
# Method 1: Using old Manifold routines.
VectorDerivative = (MF.rbasis[0]/MF.E_sq)*diff(g, u) + (MF.rbasis[1]/MF.E_sq)*diff(g, v)
print('\\eval{\\nabla g}{u=1,v=0} =', VectorDerivative.subs({u: 1, v: 0}))
# Method 2: Using new Manifold routines.
dg = MF.Grad(g)
print('\\eval{\\f{Grad}{g}}{u=1,v=0} =', dg.subs({u: 1, v: 0}))
dg = MF.grad*g
print('\\eval{\\nabla g}{u=1,v=0} =', dg.subs({u: 1, v: 0}))
return
