How to use the sympy.lambdify function in sympy
To help you get started, we’ve selected a few sympy examples, based on popular ways it is used in public projects.
Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
sympy / sympy / sympy / core / evalf.py View on Github 
quotient between successive terms must be a quotient of integer
polynomials.
"""
from sympy import Float, hypersimp, lambdify
if prec == float('inf'):
raise NotImplementedError('does not support inf prec')
if start:
expr = expr.subs(n, n + start)
hs = hypersimp(expr, n)
if hs is None:
raise NotImplementedError("a hypergeometric series is required")
num, den = hs.as_numer_denom()
func1 = lambdify(n, num)
func2 = lambdify(n, den)
h, g, p = check_convergence(num, den, n)
if h < 0:
raise ValueError("Sum diverges like (n!)^%i" % (-h))
term = expr.subs(n, 0)
if not term.is_Rational:
raise NotImplementedError("Non rational term functionality is not implemented.")
# Direct summation if geometric or faster
if h > 0 or (h == 0 and abs(g) > 1):
term = (MPZ(term.p) << prec) // term.q
s = term
k = 1r['symbolic'][0].as_symbol('C10'):15e-9,
r['symbolic'][0].as_symbol('V1'):1,
r['symbolic'][0].as_symbol('s'):1j*w,
})
out_bp = out_bp.subs({r['symbolic'][0].as_symbol('RF1'):10e3,
r['symbolic'][0].as_symbol('C10'):15e-9,
r['symbolic'][0].as_symbol('V1'):1,
r['symbolic'][0].as_symbol('s'):1j*w,
})
out_lp = out_lp.subs({r['symbolic'][0].as_symbol('RF1'):10e3,
r['symbolic'][0].as_symbol('C10'):15e-9,
r['symbolic'][0].as_symbol('V1'):1,
r['symbolic'][0].as_symbol('s'):1j*w,
})
out_lp = sympy.lambdify((w,), out_lp, modules='numpy')
out_bp = sympy.lambdify((w,), out_bp, modules='numpy')
out_hp = sympy.lambdify((w,), out_hp, modules='numpy')
ws = r['ac']['w'][::30]
fig = plt.figure()
plt.title(mycircuit.title)
plt.subplot(211)
plt.hold(True)
plt.semilogx(r['ac']['w']/2./np.pi, 20*np.log10(np.abs(r['ac']['VU1o'])), label="HP output (AC)")
plt.semilogx(r['ac']['w']/2./np.pi, 20*np.log10(np.abs(r['ac']['VU2o'])), label="BP output (AC)")
plt.semilogx(r['ac']['w']/2./np.pi, 20*np.log10(np.abs(r['ac']['VU3o'])), label="LP output (AC)")
plt.semilogx(ws/2./np.pi, 20*np.log10(np.abs(out_hp(ws))), 'v', label="HP output (SYMB)")
plt.semilogx(ws/2./np.pi, 20*np.log10(np.abs(out_bp(ws))), 'v', label="BP output (SYMB)")
plt.semilogx(ws/2./np.pi, 20*np.log10(np.abs(out_lp(ws))), 'v', label="LP output (SYMB)")
plt.hold(False)
plt.grid(True)
plt.legend()
plt.ylabel('Magnitude [dB]')from sympy import symbols, sin, lambdify
from shenfun import *
comm = MPI.COMM_WORLD
x, y = symbols("x,y")
# Some right hand side (manufactured solution)
uex = sin(2*y)*(1-y**2)
uey = sin(2*x)*(1-y**2)
pe = -0.1*sin(2*x)
fx = -uex.diff(x, 2) - uex.diff(y, 2) - pe.diff(x, 1)
fy = -uey.diff(x, 2) - uey.diff(y, 2) - pe.diff(y, 1)
h = uex.diff(x, 1) + uey.diff(y, 1)
# Lambdify for faster evaluation
ulx = lambdify((x, y), uex, 'numpy')
uly = lambdify((x, y), uey, 'numpy')
flx = lambdify((x, y), fx, 'numpy')
fly = lambdify((x, y), fy, 'numpy')
hl = lambdify((x, y), h, 'numpy')
pl = lambdify((x, y), pe, 'numpy')
N = (20, 20)
family = sys.argv[-1] if len(sys.argv) == 2 else 'Legendre'
K0 = Basis(N[0], 'Fourier', dtype='d', domain=(0, 2*np.pi))
SD = Basis(N[1], family, bc=(0, 0))
ST = Basis(N[1], family)
TD = TensorProductSpace(comm, (K0, SD), axes=(1, 0))
Q = TensorProductSpace(comm, (K0, ST), axes=(1, 0))
V = VectorTensorProductSpace(TD)
VQ = MixedTensorProductSpace([V, Q])def taylor(fx, xs, order, x_range=(0, 1), n=200):
x0, x1 = x_range
x = np.linspace(float(x0), float(x1), n)
fy = sy.lambdify(xs, fx, modules=['numpy'])(x)
tx = fx.series(xs, n=order).removeO()
if tx.is_Number:
ty = np.zeros_like(x)
ty.fill(float(tx))
else:
ty = sy.lambdify(xs, tx, modules=['numpy'])(x)
return x, fy, ty# for the n-th derivative just apply n-1th time the derivative on the result.
# However again sympy.linsolve() does not handle high dimensional case
dfds0l = sp.lambdify(symb, dfds0, modules='numpy')
dfdsil = sp.lambdify(symb, dfdsi, modules='numpy')
if di > 1:
d2fds0 = sp.diff(fexpr, symb[0], 2)
d2fdsi = sp.diff(fexpr, symb[i], 2)
d2fds0l = sp.lambdify(symb, d2fds0, modules='numpy')
d2fdsil = sp.lambdify(symb, d2fdsi, modules='numpy')
if di == 3:
d3fds0 = sp.diff(fexpr, symb[0], 3)
d3fdsi = sp.diff(fexpr, symb[i], 3)
d3fds0l = sp.lambdify(symb, d3fds0, modules='numpy')
d3fdsil = sp.lambdify(symb, d3fdsi, modules='numpy')
else:
# The implicit differentiation order can be extended if needed
error("Implicit differentiation of order n is not handled for n>3")
def multidimensional_numpy_solve(shape):
if len(shape) == 0:
return lambda x: 1/x
elif len(shape) == 2:
return np.linalg.inv
else:
return np.linalg.tensorinv
# We store the function to avoid to have a if in the loop
solve_multid = multidimensional_numpy_solve(ufl_shape)
# TODO : Vectorized version ?try:
from sympy import lambdify, symbols
except:
warn('Please install sympy to convert channel ')
X = symbols('X')
for ref in range(len(cc_ref)):
if isinstance(cc_ref[ref], CCBlock):
if cc_ref[ref]['cc_type'] == 3:
# formula to be applied
cc_ref[ref] = lambdify(X, cc_ref[ref]['cc_ref']['Comment'],
modules='numpy', dummify=False)
elif cc_ref[ref]['cc_type'] == 1: # linear conversion
cc_ref[ref] = lambdify(X, '{0} * X + {1}'.format(cc_val[1], cc_val[0]),
modules='numpy', dummify=False)
else: # identity, no conversion
cc_ref[ref] = lambdify(X, '1 * X',
modules='numpy', dummify=False)
# Otherwise a string
# look up in range keys
temp = []
for value in vector:
key_index = val_count # default index if not found
for i in range(val_count):
if key_min[i] <= value <= key_max[i]:
key_index = i
break
if callable(cc_ref[key_index]):
# TXBlock string
temp.append(cc_ref[key_index](value))
else: # scale to be applied
temp.append(cc_ref[key_index])
return asarray(temp)def eval(self, *args):
f = sympy.lambdify(self._variables, self.parametrization)
return Point(f(*args))alpha = sym_atan(rprime)
# Radius of curvature of meniscus
rm = r/sym_cos(alpha+theta)
# distance from center of curvature to meniscus contact point (Pythagoras)
d = sym_sqrt(rm**2 - r**2)
# angle between throat axis, meniscus center and meniscus contact point
gamma = sym_atan(r/d)
# Capillary Pressure
p = -2*sigma*sym_cos(alpha+theta)/r
# Callable functions
rx = lambdify((x, a, b, rt), r, 'numpy')
fill_angle = lambdify((x, a, b, rt), alpha, 'numpy')
Pc = lambdify((x, a, b, rt, sigma, theta), p, 'numpy')
rad_curve = lambdify((x, a, b, rt, sigma, theta), rm, 'numpy')
c2x = lambdify((x, a, b, rt, sigma, theta), d, 'numpy')
cap_angle = lambdify((x, a, b, rt, sigma, theta), gamma, 'numpy')
# All relative positions along throat
# pos = np.arange(-0.999, 0.999, 1/num_points)
hp = int(num_points/2)
log_pos = np.logspace(-4, -1, hp+1)[:-1]
lin_pos = np.arange(0.1, 1.0, 1/hp)
half_pos = np.concatenate((log_pos, lin_pos))
pos = np.concatenate((-half_pos[::-1], half_pos))
# Now find the positions of the menisci along each throat axis
Y, X = np.meshgrid(throatRad, pos)
X *= fa
# throat Capillary Pressure
t_Pc = Pc(X, fa, fb, Y, surface_tension, contact)
# Values of minima and maxima
Pc_min = np.min(t_Pc, axis=0)
Pc_max = np.max(t_Pc, axis=0)
# Arguments of minima and maximaeq = sp.lambdify(params,f,modules='numpy')
self.log('generated {} equations, with {} parameters'.format(
len(eqs),len(params)))
jac = None
jeqs = None
heqs = None
hessF = None
if self.algo.NeedJacobian or self.algo.NeedHessian:
# Jacobian matrix in sympy expressions
jexprs = [f.diff(x) for x in params]
if self.algo.NeedJacobian:
# Lambdified Jacobian matrix
jeqs = [sp.lambdify(params,je,modules='numpy') for je in jexprs]
self.log('generated jacobian matrix')
jac = True
if self.algo.NeedHessian:
# Lambdified Hessian matrix
heqs = [[sp.lambdify(params,je.diff(x),modules='numpy')
for x in params] for je in jexprs ]
self.log('generated hessian matrix')
hessF = self.hessF
ret = sopt.minimize(self.F,x0,(eq,jeqs,heqs), jac=jac,hess=hessF,
tol=algo.Tolerance,method=algo.getName(),options=algo.Options)
# ret = sopt.minimize(self.F,x0,(eq,None,None),method=algo.getName())
if ret.success:
for x,v in zip(params,ret.x):
param_table[x].val = vcomm = MPI.COMM_WORLD
x, y = symbols("x,y")
assert comm.Get_size() == 1, "Two non-periodic directions only have solver implemented for serial"
# Some right hand side (manufactured solution)
uex = (cos(4*np.pi*x)+sin(2*np.pi*y))*(1-y**2)*(1-x**2)
uey = (sin(2*np.pi*x)+cos(6*np.pi*y))*(1-y**2)*(1-x**2)
pe = -0.1*sin(2*x)*sin(4*y)
fx = -uex.diff(x, 2) - uex.diff(y, 2) - pe.diff(x, 1)
fy = -uey.diff(x, 2) - uey.diff(y, 2) - pe.diff(y, 1)
h = uex.diff(x, 1) + uey.diff(y, 1)
# Lambdify for faster evaluation
ulx = lambdify((x, y), uex, 'numpy')
uly = lambdify((x, y), uey, 'numpy')
flx = lambdify((x, y), fx, 'numpy')
fly = lambdify((x, y), fy, 'numpy')
hl = lambdify((x, y), h, 'numpy')
pl = lambdify((x, y), pe, 'numpy')
N = (38, 38)
#family = 'Chebyshev'
family = 'Legendre'
D0X = Basis(N[0], family, bc=(0, 0), scaled=True)
D0Y = Basis(N[1], family, bc=(0, 0), scaled=True)
PX = Basis(N[0], family)
PY = Basis(N[1], family)
# To get a P_N x P_{N-2} space, just pick the first N-2 items of the pressure basis
# Note that this effectively sets P_N and P_{N-1} to zero, but still the basis uses
# the same quadrature points as the Dirichlet basis, which is required for the inner
