def lagrange_interpolator(points, delta, delta_i_list):
"""
Construct the multidimensional Lagrange interpolating polynomial
defined by *delta* and *delta_i_list* at *points*.
This is (7) in the reference
Kamron Saniee, A Simple Expression for Multivariate Lagrange
Interpolation, SIAM, 2007.
"""
f = 0
for i, delta_i in enumerate(delta_i_list):
f += points[i][-1] * delta_i / delta
# fix for corner case when polynomial reduces to a single float
if isinstance(f, float):
return SYM.Float(f)
return f
python类Float()的实例源码
def get_paper_part_frac(degree):
from ..partialfrac import thetas_alphas_to_expr_complex
from sympy import Float, I
# Values above are negative what we expect. It also only includes one of
# each complex conjugate, which is fine so long as we pass in the positive
# imaginary parts for the thetas.
thetas = [-Float(r) + Float(i)*I for r, i in
zip(part_frac_coeffs[degree]['thetas']['real'],
part_frac_coeffs[degree]['thetas']['imaginary'])]
alphas = [-Float(r) + Float(i)*I for r, i in
zip(part_frac_coeffs[degree]['alphas']['real'],
part_frac_coeffs[degree]['alphas']['imaginary'])]
[alpha0] = [Float(r) + Float(i)*I for r, i in
zip(part_frac_coeffs[degree]['alpha0']['real'],
part_frac_coeffs[degree]['alpha0']['imaginary'])]
return thetas_alphas_to_expr_complex(thetas, alphas, alpha0)
def get_paper_thetas_alphas(degree):
from sympy import Float, I
thetas = [-Float(r) + Float(i)*I for r, i in
zip(part_frac_coeffs[degree]['thetas']['real'],
part_frac_coeffs[degree]['thetas']['imaginary'])]
thetas = thetas + [i.conjugate() for i in thetas]
alphas = [-Float(r) + Float(i)*I for r, i in
zip(part_frac_coeffs[degree]['alphas']['real'],
part_frac_coeffs[degree]['alphas']['imaginary'])]
alphas = alphas + [i.conjugate() for i in alphas]
[alpha0] = [Float(r) + Float(i)*I for r, i in
zip(part_frac_coeffs[degree]['alpha0']['real'],
part_frac_coeffs[degree]['alpha0']['imaginary'])]
return thetas, alphas, alpha0
def do_arg(a):
if isinstance(a, str):
arg = Symbol(a, integer=True)
elif isinstance(a, (Integer, Float)):
arg = a
elif isinstance(a, ArithmeticExpression):
arg = a.expr
if isinstance(arg, (Symbol, Variable)):
arg = Symbol(arg.name, integer=True)
else:
arg = convert_to_integer_expression(arg)
else:
raise Exception('Wrong instance in do_arg')
return arg
# ...
# ... TODO improve. this version is not working with function calls
def extract_arg(self, name):
"""
returns an argument as a variable, given its name
name: str
variable name
"""
if name is None:
return None
var = None
if isinstance(name, (Integer, Float)):
var = Integer(name)
elif isinstance(name, str):
if name in namespace:
var = namespace[name]
else:
raise Exception("could not find {} in namespace ".format(name))
elif isinstance(name, ArithmeticExpression):
var = do_arg(name)
else:
raise Exception("Unexpected type {0} for {1}".format(type(name), name))
return var
def test_precision():
f1 = Float("1.01")
f2 = Float("1.0000000000000000000001")
for x in [1, 1e-2, 1e-6, 1e-13, 1e-14, 1e-16, 1e-20, 1e-100, 1e-300,
f1, f2]:
result = construct_domain([x])
y = float(result[1][0])
assert abs(x - y) / x < 1e-14 # Test relative accuracy
result = construct_domain([f1])
y = result[1][0]
assert y-1 > 1e-50
result = construct_domain([f2])
y = result[1][0]
assert y-1 > 1e-50
def test_pow():
n_x = f_x**3
assert n_x.ndim == 1
assert isinstance(n_x, Formula1D)
assert_allclose(n_x(x).v, (f_x(x)**3).v)
assert str(n_x(x).units) == "km**3"
m_x = f_x**0.5
assert n_x.ndim == 1
assert isinstance(m_x, Formula1D)
assert_allclose(m_x(x).v, (f_x(x)**0.5).v)
assert str(m_x(x).units) == "sqrt(km)"
l_x = f_x**Float(0.3)
assert l_x.ndim == 1
assert isinstance(l_x, Formula1D)
assert_allclose(l_x(x).v, (f_x(x)**0.3).v)
assert str(l_x(x).units) == "km**(3/10)"
def _gauss_from_coefficients_mpmath(alpha, beta, decimal_places):
mp.dps = decimal_places
# Create vector cut of the first value of beta
n = len(alpha)
b = mp.zeros(n, 1)
for i in range(n-1):
b[i] = mp.sqrt(beta[i+1])
z = mp.zeros(1, n)
z[0, 0] = 1
d = mp.matrix(alpha)
tridiag_eigen(mp, d, b, z)
# nx1 matrix -> list of sympy floats
x = numpy.array([sympy.Float(xx) for xx in d])
w = numpy.array([beta[0] * mp.power(ww, 2) for ww in z])
return x, w
def format_kinetics(self, rate = False, precision = 3):
"""Convert the kinetic parameter or rate to string.
If rate = True, return a string representing the rate,
otherwise one representing the kinetic parameter.
If the kinetic parameter is a float, use exponent notation,
with a number of digits equal to precision (3 is the default).
If the reaction has no defined rate, return None."""
k = None
if self.rate:
if isinstance(self.kinetic_param, sp.Float):
k = "{:.{}e}".format(self.kinetic_param, precision)
if rate:
k = k + "*" + str(self.reactant.ma())
else:
if rate:
k = str(self.rate)
else:
k = str(self.kinetic_param)
return k
def create_expression(n):
if n not in coeffs:
raise ValueError("Don't have coefficients for {}".format(n))
from sympy import Float, Add
p = coeffs[n]['p']
q = coeffs[n]['q']
# Don't use Poly here, it loses precision
num = Add(*[Float(p[i])*t**i for i in range(n+1)])
den = Add(*[Float(q[i])*t**i for i in range(n+1)])
return num/den
def test_scalar_sympy():
assert represent(Integer(1)) == Integer(1)
assert represent(Float(1.0)) == Float(1.0)
assert represent(1.0 + I) == 1.0 + I
def test_scalar_numpy():
if not np:
skip("numpy not installed.")
assert represent(Integer(1), format='numpy') == 1
assert represent(Float(1.0), format='numpy') == 1.0
assert represent(1.0 + I, format='numpy') == 1.0 + 1.0j
def test_scalar_scipy_sparse():
if not np:
skip("numpy not installed.")
if not scipy:
skip("scipy not installed.")
assert represent(Integer(1), format='scipy.sparse') == 1
assert represent(Float(1.0), format='scipy.sparse') == 1.0
assert represent(1.0 + I, format='scipy.sparse') == 1.0 + 1.0j
def test_hbar():
assert hbar.is_commutative is True
assert hbar.is_real is True
assert hbar.is_positive is True
assert hbar.is_negative is False
assert hbar.is_irrational is True
assert hbar.evalf() == Float(1.05457162e-34)
def __new__(cls, variable, value):
if not isinstance(variable, Variable):
raise TypeError("variable must be of type Variable")
_valid_instances = (Nil, Variable,
IndexedVariable, IndexedElement,
Tuple,
int, float, bool, complex, str,
Boolean, sp_Integer, sp_Float)
if not isinstance(value, _valid_instances):
raise TypeError('non-valid instance for value, '
'given {0}'.format(type(value)))
return Basic.__new__(cls, variable, value)
def is_simple_assign(expr):
if not isinstance(expr, Assign):
return False
assignable = [Variable, IndexedVariable, IndexedElement]
assignable += [sp_Integer, sp_Float]
assignable = tuple(assignable)
if isinstance(expr.rhs, assignable):
return True
else:
return False
def test_scalar_sympy():
assert represent(Integer(1)) == Integer(1)
assert represent(Float(1.0)) == Float(1.0)
assert represent(1.0 + I) == 1.0 + I
def test_scalar_numpy():
if not np:
skip("numpy not installed.")
assert represent(Integer(1), format='numpy') == 1
assert represent(Float(1.0), format='numpy') == 1.0
assert represent(1.0 + I, format='numpy') == 1.0 + 1.0j
def test_scalar_scipy_sparse():
if not np:
skip("numpy not installed.")
if not scipy:
skip("scipy not installed.")
assert represent(Integer(1), format='scipy.sparse') == 1
assert represent(Float(1.0), format='scipy.sparse') == 1.0
assert represent(1.0 + I, format='scipy.sparse') == 1.0 + 1.0j
def test_hbar():
assert hbar.is_commutative is True
assert hbar.is_real is True
assert hbar.is_positive is True
assert hbar.is_negative is False
assert hbar.is_irrational is True
assert hbar.evalf() == Float(1.05457162e-34)
def feq(a, b):
"""Test if two floating point values are 'equal'."""
t_float = Float("1.0E-10")
return -t_float < a - b < t_float
def feq(a, b):
"""Test if two floating point values are 'equal'."""
t_float = Float("1.0E-10")
return -t_float < a - b < t_float
def E_nl_dirac(n, l, spin_up=True, Z=1, c=Float("137.035999037")):
"""
Returns the relativistic energy of the state (n, l, spin) in Hartree atomic
units.
The energy is calculated from the Dirac equation. The rest mass energy is
*not* included.
n, l
quantum numbers 'n' and 'l'
spin_up
True if the electron spin is up (default), otherwise down
Z
atomic number (1 for Hydrogen, 2 for Helium, ...)
c
speed of light in atomic units. Default value is 137.035999037,
taken from: http://arxiv.org/abs/1012.3627
Examples
========
>>> from sympy.physics.hydrogen import E_nl_dirac
>>> E_nl_dirac(1, 0)
-0.500006656595360
>>> E_nl_dirac(2, 0)
-0.125002080189006
>>> E_nl_dirac(2, 1)
-0.125000416028342
>>> E_nl_dirac(2, 1, False)
-0.125002080189006
>>> E_nl_dirac(3, 0)
-0.0555562951740285
>>> E_nl_dirac(3, 1)
-0.0555558020932949
>>> E_nl_dirac(3, 1, False)
-0.0555562951740285
>>> E_nl_dirac(3, 2)
-0.0555556377366884
>>> E_nl_dirac(3, 2, False)
-0.0555558020932949
"""
if not (l >= 0):
raise ValueError("'l' must be positive or zero")
if not (n > l):
raise ValueError("'n' must be greater than 'l'")
if (l == 0 and spin_up is False):
raise ValueError("Spin must be up for l==0.")
# skappa is sign*kappa, where sign contains the correct sign
if spin_up:
skappa = -l - 1
else:
skappa = -l
c = S(c)
beta = sqrt(skappa**2 - Z**2/c**2)
return c**2/sqrt(1 + Z**2/(n + skappa + beta)**2/c**2) - c**2
def E_nl_dirac(n, l, spin_up=True, Z=1, c=Float("137.035999037")):
"""
Returns the relativistic energy of the state (n, l, spin) in Hartree atomic
units.
The energy is calculated from the Dirac equation. The rest mass energy is
*not* included.
n, l
quantum numbers 'n' and 'l'
spin_up
True if the electron spin is up (default), otherwise down
Z
atomic number (1 for Hydrogen, 2 for Helium, ...)
c
speed of light in atomic units. Default value is 137.035999037,
taken from: http://arxiv.org/abs/1012.3627
Examples
========
>>> from sympy.physics.hydrogen import E_nl_dirac
>>> E_nl_dirac(1, 0)
-0.500006656595360
>>> E_nl_dirac(2, 0)
-0.125002080189006
>>> E_nl_dirac(2, 1)
-0.125000416028342
>>> E_nl_dirac(2, 1, False)
-0.125002080189006
>>> E_nl_dirac(3, 0)
-0.0555562951740285
>>> E_nl_dirac(3, 1)
-0.0555558020932949
>>> E_nl_dirac(3, 1, False)
-0.0555562951740285
>>> E_nl_dirac(3, 2)
-0.0555556377366884
>>> E_nl_dirac(3, 2, False)
-0.0555558020932949
"""
if not (l >= 0):
raise ValueError("'l' must be positive or zero")
if not (n > l):
raise ValueError("'n' must be greater than 'l'")
if (l == 0 and spin_up is False):
raise ValueError("Spin must be up for l==0.")
# skappa is sign*kappa, where sign contains the correct sign
if spin_up:
skappa = -l - 1
else:
skappa = -l
c = S(c)
beta = sqrt(skappa**2 - Z**2/c**2)
return c**2/sqrt(1 + Z**2/(n + skappa + beta)**2/c**2) - c**2
