我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用sympy.cos()。
def __init__(self, n): self.degree = 2 self.dim = n n2 = fr(n, 2) if n % 2 == 0 else fr(n-1, 2) pts = [[ [sqrt(fr(2, n+2)) * cos(2*k*i*pi / (n+1)) for i in range(n+1)], [sqrt(fr(2, n+2)) * sin(2*k*i*pi / (n+1)) for i in range(n+1)], ] for k in range(1, n2+1)] if n % 2 == 1: sqrt3pm = numpy.full(n+1, 1/sqrt(n+2)) sqrt3pm[1::2] *= -1 pts.append(sqrt3pm) pts = numpy.vstack(pts).T data = [(fr(1, n+1), pts)] self.points, self.weights = untangle(data) self.weights *= volume_unit_ball(n) return
def test_lie_derivative(self): Lfh = pm.lie_derivatives(self.h, self.f, self.x, 0) self.assertEqual(Lfh, [self.h]) Lfh = pm.lie_derivatives(self.h, self.f, self.x, 1) self.assertEqual(Lfh, [self.h, sp.Matrix([-2*self._x1*self._x2**2 - sp.sin(self._x1)*sp.cos(self._x2)]) ]) Lfh = pm.lie_derivatives(self.h, self.f, self.x, 2) self.assertEqual(Lfh, [self.h, sp.Matrix([-2*self._x1*self._x2**2 - sp.sin(self._x1)*sp.cos(self._x2)]), sp.Matrix([ -self._x2**2*( -2*self._x2**2 - sp.cos(self._x1)*sp.cos(self._x2) ) + sp.sin(self._x1)*( - 4*self._x1*self._x2 + sp.sin(self._x1)*sp.sin(self._x2) )]) ])
def test_epath_apply(): expr = [((x, 1, t), 2), ((3, y, 4), z)] func = lambda expr: expr**2 assert epath("/*", expr, list) == [[(x, 1, t), 2], [(3, y, 4), z]] assert epath("/*/[0]", expr, list) == [([x, 1, t], 2), ([3, y, 4], z)] assert epath("/*/[1]", expr, func) == [((x, 1, t), 4), ((3, y, 4), z**2)] assert epath("/*/[2]", expr, list) == expr assert epath("/*/[0]/int", expr, func) == [((x, 1, t), 2), ((9, y, 16), z)] assert epath("/*/[0]/Symbol", expr, func) == [((x**2, 1, t**2), 2), ((3, y**2, 4), z)] assert epath( "/*/[0]/int[1:]", expr, func) == [((x, 1, t), 2), ((3, y, 16), z)] assert epath("/*/[0]/Symbol[1:]", expr, func) == [((x, 1, t**2), 2), ((3, y**2, 4), z)] assert epath("/Symbol", x + y + z + 1, func) == x**2 + y**2 + z**2 + 1 assert epath("/*/*/Symbol", t + sin(x + 1) + cos(x + y + E), func) == \ t + sin(x**2 + 1) + cos(x**2 + y**2 + E)
def test_sqrt_symbolic_denest(): x = Symbol('x') z = sqrt(((1 + sqrt(sqrt(2 + x) + 3))**2).expand()) assert sqrtdenest(z) == sqrt((1 + sqrt(sqrt(2 + x) + 3))**2) z = sqrt(((1 + sqrt(sqrt(2 + cos(1)) + 3))**2).expand()) assert sqrtdenest(z) == 1 + sqrt(sqrt(2 + cos(1)) + 3) z = ((1 + cos(2))**4 + 1).expand() assert sqrtdenest(z) == z z = sqrt(((1 + sqrt(sqrt(2 + cos(3*x)) + 3))**2 + 1).expand()) assert sqrtdenest(z) == z c = cos(3) c2 = c**2 assert sqrtdenest(sqrt(2*sqrt(1 + r3)*c + c2 + 1 + r3*c2)) == \ -1 - sqrt(1 + r3)*c ra = sqrt(1 + r3) z = sqrt(20*ra*sqrt(3 + 3*r3) + 12*r3*ra*sqrt(3 + 3*r3) + 64*r3 + 112) assert sqrtdenest(z) == z
def equation(self, type='cosine'): """ Returns equation of the wave. Examples ======== >>> from sympy import symbols >>> from sympy.physics.optics import TWave >>> A, phi, f = symbols('A, phi, f') >>> w = TWave(A, f, phi) >>> w.equation('cosine') A*cos(2*pi*f*t + phi) """ if not isinstance(type, str): raise TypeError("type can only be a string.") if type == 'cosine': return self._amplitude*cos(self.angular_velocity*Symbol('t') + self._phase)
def test_twave(): A1, phi1, A2, phi2, f = symbols('A1, phi1, A2, phi2, f') c = Symbol('c') # Speed of light in vacuum n = Symbol('n') # Refractive index t = Symbol('t') # Time w1 = TWave(A1, f, phi1) w2 = TWave(A2, f, phi2) assert w1.amplitude == A1 assert w1.frequency == f assert w1.phase == phi1 assert w1.wavelength == c/(f*n) assert w1.time_period == 1/f w3 = w1 + w2 assert w3.amplitude == sqrt(A1**2 + 2*A1*A2*cos(phi1 - phi2) + A2**2) assert w3.frequency == f assert w3.wavelength == c/(f*n) assert w3.time_period == 1/f assert w3.angular_velocity == 2*pi*f assert w3.equation() == sqrt(A1**2 + 2*A1*A2*cos(phi1 - phi2) + A2**2)*cos(2*pi*f*t + phi1 + phi2)
def test_pend(): q, u = dynamicsymbols('q u') qd, ud = dynamicsymbols('q u', 1) m, l, g = symbols('m l g') N = ReferenceFrame('N') P = Point('P') P.set_vel(N, -l * u * sin(q) * N.x + l * u * cos(q) * N.y) kd = [qd - u] FL = [(P, m * g * N.x)] pa = Particle('pa', P, m) BL = [pa] KM = KanesMethod(N, [q], [u], kd) KM.kanes_equations(FL, BL) MM = KM.mass_matrix forcing = KM.forcing rhs = MM.inv() * forcing rhs.simplify() assert expand(rhs[0]) == expand(-g / l * sin(q))
def test_dot_different_frames(): assert dot(N.x, A.x) == cos(q1) assert dot(N.x, A.y) == -sin(q1) assert dot(N.x, A.z) == 0 assert dot(N.y, A.x) == sin(q1) assert dot(N.y, A.y) == cos(q1) assert dot(N.y, A.z) == 0 assert dot(N.z, A.x) == 0 assert dot(N.z, A.y) == 0 assert dot(N.z, A.z) == 1 assert dot(N.x, A.x + A.y) == sqrt(2)*cos(q1 + pi/4) == dot(A.x + A.y, N.x) assert dot(A.x, C.x) == cos(q3) assert dot(A.x, C.y) == 0 assert dot(A.x, C.z) == sin(q3) assert dot(A.y, C.x) == sin(q2)*sin(q3) assert dot(A.y, C.y) == cos(q2) assert dot(A.y, C.z) == -sin(q2)*cos(q3) assert dot(A.z, C.x) == -cos(q2)*sin(q3) assert dot(A.z, C.y) == sin(q2) assert dot(A.z, C.z) == cos(q2)*cos(q3)
def test_cross_different_frames(): assert cross(N.x, A.x) == sin(q1)*A.z assert cross(N.x, A.y) == cos(q1)*A.z assert cross(N.x, A.z) == -sin(q1)*A.x - cos(q1)*A.y assert cross(N.y, A.x) == -cos(q1)*A.z assert cross(N.y, A.y) == sin(q1)*A.z assert cross(N.y, A.z) == cos(q1)*A.x - sin(q1)*A.y assert cross(N.z, A.x) == A.y assert cross(N.z, A.y) == -A.x assert cross(N.z, A.z) == 0 assert cross(N.x, A.x) == sin(q1)*A.z assert cross(N.x, A.y) == cos(q1)*A.z assert cross(N.x, A.x + A.y) == sin(q1)*A.z + cos(q1)*A.z assert cross(A.x + A.y, N.x) == -sin(q1)*A.z - cos(q1)*A.z assert cross(A.x, C.x) == sin(q3)*C.y assert cross(A.x, C.y) == -sin(q3)*C.x + cos(q3)*C.z assert cross(A.x, C.z) == -cos(q3)*C.y assert cross(C.x, A.x) == -sin(q3)*C.y assert cross(C.y, A.x) == sin(q3)*C.x - cos(q3)*C.z assert cross(C.z, A.x) == cos(q3)*C.y
def test_partial_velocity(): q1, q2, q3, u1, u2, u3 = dynamicsymbols('q1 q2 q3 u1 u2 u3') u4, u5 = dynamicsymbols('u4, u5') r = symbols('r') N = ReferenceFrame('N') Y = N.orientnew('Y', 'Axis', [q1, N.z]) L = Y.orientnew('L', 'Axis', [q2, Y.x]) R = L.orientnew('R', 'Axis', [q3, L.y]) R.set_ang_vel(N, u1 * L.x + u2 * L.y + u3 * L.z) C = Point('C') C.set_vel(N, u4 * L.x + u5 * (Y.z ^ L.x)) Dmc = C.locatenew('Dmc', r * L.z) Dmc.v2pt_theory(C, N, R) vel_list = [Dmc.vel(N), C.vel(N), R.ang_vel_in(N)] u_list = [u1, u2, u3, u4, u5] assert (partial_velocity(vel_list, u_list, N) == [[- r*L.y, r*L.x, 0, L.x, cos(q2)*L.y - sin(q2)*L.z], [0, 0, 0, L.x, cos(q2)*L.y - sin(q2)*L.z], [L.x, L.y, L.z, 0, 0]])
def test_expand(): m0 = Matrix([[x*(x + y), 2], [((x + y)*y)*x, x*(y + x*(x + y))]]) # Test if expand() returns a matrix m1 = m0.expand() assert m1 == Matrix( [[x*y + x**2, 2], [x*y**2 + y*x**2, x*y + y*x**2 + x**3]]) a = Symbol('a', real=True) assert Matrix([exp(I*a)]).expand(complex=True) == \ Matrix([cos(a) + I*sin(a)]) assert Matrix([[0, 1, 2], [0, 0, -1], [0, 0, 0]]).exp() == Matrix([ [1, 1, Rational(3, 2)], [0, 1, -1], [0, 0, 1]] )
def test_simplify(): f, n = symbols('f, n') m = Matrix([[1, x], [x + 1/x, x - 1]]) m = m.row_join(eye(m.cols)) raw = m.rref(simplify=lambda x: x)[0] assert raw != m.rref(simplify=True)[0] M = Matrix([[ 1/x + 1/y, (x + x*y) / x ], [ (f(x) + y*f(x))/f(x), 2 * (1/n - cos(n * pi)/n) / pi ]]) M.simplify() assert M == Matrix([[ (x + y)/(x * y), 1 + y ], [ 1 + y, 2*((1 - 1*cos(pi*n))/(pi*n)) ]]) eq = (1 + x)**2 M = Matrix([[eq]]) M.simplify() assert M == Matrix([[eq]]) M.simplify(ratio=oo) == M assert M == Matrix([[eq.simplify(ratio=oo)]])
def test_rotation_matrices(): # This tests the rotation matrices by rotating about an axis and back. theta = pi/3 r3_plus = rot_axis3(theta) r3_minus = rot_axis3(-theta) r2_plus = rot_axis2(theta) r2_minus = rot_axis2(-theta) r1_plus = rot_axis1(theta) r1_minus = rot_axis1(-theta) assert r3_minus*r3_plus*eye(3) == eye(3) assert r2_minus*r2_plus*eye(3) == eye(3) assert r1_minus*r1_plus*eye(3) == eye(3) # Check the correctness of the trace of the rotation matrix assert r1_plus.trace() == 1 + 2*cos(theta) assert r2_plus.trace() == 1 + 2*cos(theta) assert r3_plus.trace() == 1 + 2*cos(theta) # Check that a rotation with zero angle doesn't change anything. assert rot_axis1(0) == eye(3) assert rot_axis2(0) == eye(3) assert rot_axis3(0) == eye(3)
def test_Znm(): # http://en.wikipedia.org/wiki/Solid_harmonics#List_of_lowest_functions th, ph = Symbol("theta", real=True), Symbol("phi", real=True) from sympy.abc import n,m assert Znm(0, 0, th, ph) == Ynm(0, 0, th, ph) assert Znm(1, -1, th, ph) == (-sqrt(2)*I*(Ynm(1, 1, th, ph) - exp(-2*I*ph)*Ynm(1, 1, th, ph))/2) assert Znm(1, 0, th, ph) == Ynm(1, 0, th, ph) assert Znm(1, 1, th, ph) == (sqrt(2)*(Ynm(1, 1, th, ph) + exp(-2*I*ph)*Ynm(1, 1, th, ph))/2) assert Znm(0, 0, th, ph).expand(func=True) == 1/(2*sqrt(pi)) assert Znm(1, -1, th, ph).expand(func=True) == (sqrt(3)*I*sqrt(-cos(th)**2 + 1)*exp(I*ph)/(4*sqrt(pi)) - sqrt(3)*I*sqrt(-cos(th)**2 + 1)*exp(-I*ph)/(4*sqrt(pi))) assert Znm(1, 0, th, ph).expand(func=True) == sqrt(3)*cos(th)/(2*sqrt(pi)) assert Znm(1, 1, th, ph).expand(func=True) == (-sqrt(3)*sqrt(-cos(th)**2 + 1)*exp(I*ph)/(4*sqrt(pi)) - sqrt(3)*sqrt(-cos(th)**2 + 1)*exp(-I*ph)/(4*sqrt(pi))) assert Znm(2, -1, th, ph).expand(func=True) == (sqrt(15)*I*sqrt(-cos(th)**2 + 1)*exp(I*ph)*cos(th)/(4*sqrt(pi)) - sqrt(15)*I*sqrt(-cos(th)**2 + 1)*exp(-I*ph)*cos(th)/(4*sqrt(pi))) assert Znm(2, 0, th, ph).expand(func=True) == 3*sqrt(5)*cos(th)**2/(4*sqrt(pi)) - sqrt(5)/(4*sqrt(pi)) assert Znm(2, 1, th, ph).expand(func=True) == (-sqrt(15)*sqrt(-cos(th)**2 + 1)*exp(I*ph)*cos(th)/(4*sqrt(pi)) - sqrt(15)*sqrt(-cos(th)**2 + 1)*exp(-I*ph)*cos(th)/(4*sqrt(pi)))
def test_function_series1(): """Create our new "sin" function.""" class my_function(Function): def fdiff(self, argindex=1): return cos(self.args[0]) @classmethod def eval(cls, arg): arg = sympify(arg) if arg == 0: return sympify(0) #Test that the taylor series is correct assert my_function(x).series(x, 0, 10) == sin(x).series(x, 0, 10) assert limit(my_function(x)/x, x, 0) == 1
def test_function_series2(): """Create our new "cos" function.""" class my_function2(Function): def fdiff(self, argindex=1): return -sin(self.args[0]) @classmethod def eval(cls, arg): arg = sympify(arg) if arg == 0: return sympify(1) #Test that the taylor series is correct assert my_function2(x).series(x, 0, 10) == cos(x).series(x, 0, 10)
def test_piecewise_fold_piecewise_in_cond(): p1 = Piecewise((cos(x), x < 0), (0, True)) p2 = Piecewise((0, Eq(p1, 0)), (p1 / Abs(p1), True)) p3 = piecewise_fold(p2) assert(p2.subs(x, -pi/2) == 0.0) assert(p2.subs(x, 1) == 0.0) assert(p2.subs(x, -pi/4) == 1.0) p4 = Piecewise((0, Eq(p1, 0)), (1,True)) assert(piecewise_fold(p4) == Piecewise( (0, Or(And(Eq(cos(x), 0), x < 0), Not(x < 0))), (1, True))) r1 = 1 < Piecewise((1, x < 1), (3, True)) assert(piecewise_fold(r1) == Not(x < 1)) p5 = Piecewise((1, x < 0), (3, True)) p6 = Piecewise((1, x < 1), (3, True)) p7 = piecewise_fold(Piecewise((1, p5 < p6), (0, True))) assert(Piecewise((1, And(Not(x < 1), x < 0)), (0, True)))
def _compute_transform(self, F, s, x, **hints): from sympy import postorder_traversal global _allowed if _allowed is None: from sympy import ( exp, gamma, sin, cos, tan, cot, cosh, sinh, tanh, coth, factorial, rf) _allowed = set( [exp, gamma, sin, cos, tan, cot, cosh, sinh, tanh, coth, factorial, rf]) for f in postorder_traversal(F): if f.is_Function and f.has(s) and f.func not in _allowed: raise IntegralTransformError('Inverse Mellin', F, 'Component %s not recognised.' % f) strip = self.fundamental_strip return _inverse_mellin_transform(F, s, x, strip, **hints)
def _sine_cosine_transform(f, x, k, a, b, K, name, simplify=True): """ Compute a general sine or cosine-type transform F(k) = a int_0^oo b*sin(x*k) f(x) dx. F(k) = a int_0^oo b*cos(x*k) f(x) dx. For suitable choice of a and b, this reduces to the standard sine/cosine and inverse sine/cosine transforms. """ F = integrate(a*f*K(b*x*k), (x, 0, oo)) if not F.has(Integral): return _simplify(F, simplify), True if not F.is_Piecewise: raise IntegralTransformError(name, f, 'could not compute integral') F, cond = F.args[0] if F.has(Integral): raise IntegralTransformError(name, f, 'integral in unexpected form') return _simplify(F, simplify), cond
def trig_tansec_rule(integral): integrand, symbol = integral integrand = integrand.subs({ 1 / sympy.cos(symbol): sympy.sec(symbol) }) if any(integrand.has(f) for f in (sympy.tan, sympy.sec)): pattern, a, b, m, n = tansec_pattern(symbol) match = integrand.match(pattern) if match: a, b, m, n = match.get(a, 0),match.get(b, 0), match.get(m, 0), match.get(n, 0) return multiplexer({ tansec_tanodd_condition: tansec_tanodd, tansec_seceven_condition: tansec_seceven, tan_tansquared_condition: tan_tansquared })((a, b, m, n, integrand, symbol))
def trig_cotcsc_rule(integral): integrand, symbol = integral integrand = integrand.subs({ 1 / sympy.sin(symbol): sympy.csc(symbol), 1 / sympy.tan(symbol): sympy.cot(symbol), sympy.cos(symbol) / sympy.tan(symbol): sympy.cot(symbol) }) if any(integrand.has(f) for f in (sympy.cot, sympy.csc)): pattern, a, b, m, n = cotcsc_pattern(symbol) match = integrand.match(pattern) if match: a, b, m, n = match.get(a, 0),match.get(b, 0), match.get(m, 0), match.get(n, 0) return multiplexer({ cotcsc_cotodd_condition: cotcsc_cotodd, cotcsc_csceven_condition: cotcsc_csceven })((a, b, m, n, integrand, symbol))
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_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 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 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 main(): x = sympy.Symbol('x') y = sympy.Symbol('y') e = 1/sympy.cos(x) print() pprint(e) print('\n') pprint(e.subs(sympy.cos(x), y)) print('\n') pprint(e.subs(sympy.cos(x), y).subs(y, x**2)) e = 1/sympy.log(x) e = e.subs(x, sympy.Float("2.71828")) print('\n') pprint(e) print('\n') pprint(e.evalf()) print()
def real(self, nested_scope=None): """Return the correspond floating point number.""" op = self.children[0].name expr = self.children[1] dispatch = { 'sin': sympy.sin, 'cos': sympy.cos, 'tan': sympy.tan, 'asin': sympy.asin, 'acos': sympy.acos, 'atan': sympy.atan, 'exp': sympy.exp, 'ln': sympy.log, 'sqrt': sympy.sqrt } if op in dispatch: arg = expr.real(nested_scope) return dispatch[op](arg) else: raise NodeException("internal error: undefined external")
def sym(self, nested_scope=None): """Return the corresponding symbolic expression.""" op = self.children[0].name expr = self.children[1] dispatch = { 'sin': sympy.sin, 'cos': sympy.cos, 'tan': sympy.tan, 'asin': sympy.asin, 'acos': sympy.acos, 'atan': sympy.atan, 'exp': sympy.exp, 'ln': sympy.log, 'sqrt': sympy.sqrt } if op in dispatch: arg = expr.sym(nested_scope) return dispatch[op](arg) else: raise NodeException("internal error: undefined external")
def two_wheel_2d_model(x, u, dt): """Two wheel 2D motion model Parameters ---------- x : u : dt : Returns ------- """ gdot = np.array([[u[0, 0] * cos(x[2, 0]) * dt], [u[0, 0] * sin(x[2, 0]) * dt], [u[1, 0] * dt]]) return x + gdot
def two_wheel_3d_model(x, u, dt): """Two wheel 3D motion model Parameters ---------- x : u : dt : Returns ------- """ g1 = x[0, 0] + u[0] * cos(x[3, 0]) * dt g2 = x[1, 0] + u[0] * sin(x[3, 0]) * dt g3 = x[2, 0] + u[1] * dt g4 = x[3, 0] + u[2] * dt return np.array([g1, g2, g3, g4])
def two_wheel_3d_deriv(): """ """ x1, x2, x3, x4, x5, x6, x7 = sympy.symbols("x1,x2,x3,x4,x5,x6,x7") dt = sympy.symbols("dt") # x1 - x # x2 - y # x3 - z # x4 - theta # x5 - v # x6 - omega # x7 - vz # x, y, z, theta, v, omega, vz f1 = x1 + x5 * sympy.cos(x4) * dt f2 = x2 + x5 * sympy.sin(x4) * dt f3 = x3 + x7 * dt f4 = x4 + x6 * dt f5 = x5 f6 = x6 f7 = x7 F = sympy.Matrix([f1, f2, f3, f4, f5, f6, f7]) pprint(F.jacobian([x1, x2, x3, x4, x5, x6, x7]))
def linearity_index_inverse_depth(): """Linearity index of Inverse Depth Parameterization""" D, rho, rho_0, d_1, sigma_rho = sympy.symbols("D,rho,rho_0,d_1,sigma_rho") alpha = sympy.symbols("alpha") u = (rho * sympy.sin(alpha)) / (rho_0 * d_1 * (rho_0 - rho) + rho * sympy.cos(alpha)) # NOQA # first order derivative of u u_p = sympy.diff(u, rho) u_p = sympy.simplify(u_p) # second order derivative of u u_pp = sympy.diff(u_p, rho) u_pp = sympy.simplify(u_pp) # Linearity index L = (u_pp * 2 * sigma_rho) / (u_p) L = sympy.simplify(L) print() print("u: ", u) print("u': ", u_p) print("u'': ", u_pp) # print("L = ", L) print("L = ", L.subs(rho, 0)) print()
def compute_stencil(order, n, x_value, h_value, x0=0.): """ computes a stencil of Order order """ h,x = symbols('h x') xs = [x+h*cos(i*pi/n) for i in range(n,-1,-1)] cs = finite_diff_weights(order, xs, x0)[order][n] cs = [simplify(c) for c in cs] cs = [simplify(expr.subs(h, h_value)) for expr in cs] xs = [simplify(expr.subs(h, h_value)) for expr in xs] cs = [simplify(expr.subs(x, x_value)) for expr in cs] xs = [simplify(expr.subs(x, x_value)) for expr in xs] return xs, cs ##############################################
def test_pow_eval(): # XXX Pow does not fully support conversion of negative numbers # to their complex equivalent assert sqrt(-1) == I assert sqrt(-4) == 2*I assert sqrt( 4) == 2 assert (8)**Rational(1, 3) == 2 assert (-8)**Rational(1, 3) == 2*((-1)**Rational(1, 3)) assert sqrt(-2) == I*sqrt(2) assert (-1)**Rational(1, 3) != I assert (-10)**Rational(1, 3) != I*((10)**Rational(1, 3)) assert (-2)**Rational(1, 4) != (2)**Rational(1, 4) assert 64**Rational(1, 3) == 4 assert 64**Rational(2, 3) == 16 assert 24/sqrt(64) == 3 assert (-27)**Rational(1, 3) == 3*(-1)**Rational(1, 3) assert (cos(2) / tan(2))**2 == (cos(2) / tan(2))**2
def _eval_power(self, expt): """ b is I = sqrt(-1) e is symbolic object but not equal to 0, 1 I**r -> (-1)**(r/2) -> exp(r/2*Pi*I) -> sin(Pi*r/2) + cos(Pi*r/2)*I, r is decimal I**0 mod 4 -> 1 I**1 mod 4 -> I I**2 mod 4 -> -1 I**3 mod 4 -> -I """ if isinstance(expt, Number): if isinstance(expt, Integer): expt = expt.p % 4 if expt == 0: return S.One if expt == 1: return S.ImaginaryUnit if expt == 2: return -S.One return -S.ImaginaryUnit return (S.NegativeOne)**(expt*S.Half) return
def test_condition_number(): x = Symbol('x', real=True) A = eye(3) A[0, 0] = 10 A[2, 2] = S(1)/10 assert A.condition_number() == 100 A[1, 1] = x assert A.condition_number() == Max(10, Abs(x)) / Min(S(1)/10, Abs(x)) M = Matrix([[cos(x), sin(x)], [-sin(x), cos(x)]]) Mc = M.condition_number() assert all(Float(1.).epsilon_eq(Mc.subs(x, val).evalf()) for val in [Rational(1, 5), Rational(1, 2), Rational(1, 10), pi/2, pi, 7*pi/4 ]) #issue 10782 assert Matrix([]).condition_number() == 0
def test_derivative_by_array(): from sympy.abc import a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z bexpr = x*y**2*exp(z)*log(t) sexpr = sin(bexpr) cexpr = cos(bexpr) a = Array([sexpr]) assert derive_by_array(sexpr, t) == x*y**2*exp(z)*cos(x*y**2*exp(z)*log(t))/t assert derive_by_array(sexpr, [x, y, z]) == Array([bexpr/x*cexpr, 2*y*bexpr/y**2*cexpr, bexpr*cexpr]) assert derive_by_array(a, [x, y, z]) == Array([[bexpr/x*cexpr], [2*y*bexpr/y**2*cexpr], [bexpr*cexpr]]) assert derive_by_array(sexpr, [[x, y], [z, t]]) == Array([[bexpr/x*cexpr, 2*y*bexpr/y**2*cexpr], [bexpr*cexpr, bexpr/log(t)/t*cexpr]]) assert derive_by_array(a, [[x, y], [z, t]]) == Array([[[bexpr/x*cexpr], [2*y*bexpr/y**2*cexpr]], [[bexpr*cexpr], [bexpr/log(t)/t*cexpr]]]) assert derive_by_array([[x, y], [z, t]], [x, y]) == Array([[[1, 0], [0, 0]], [[0, 1], [0, 0]]]) assert derive_by_array([[x, y], [z, t]], [[x, y], [z, t]]) == Array([[[[1, 0], [0, 0]], [[0, 1], [0, 0]]], [[[0, 0], [1, 0]], [[0, 0], [0, 1]]]])
def test_Znm(): # http://en.wikipedia.org/wiki/Solid_harmonics#List_of_lowest_functions th, ph = Symbol("theta", real=True), Symbol("phi", real=True) assert Znm(0, 0, th, ph) == Ynm(0, 0, th, ph) assert Znm(1, -1, th, ph) == (-sqrt(2)*I*(Ynm(1, 1, th, ph) - exp(-2*I*ph)*Ynm(1, 1, th, ph))/2) assert Znm(1, 0, th, ph) == Ynm(1, 0, th, ph) assert Znm(1, 1, th, ph) == (sqrt(2)*(Ynm(1, 1, th, ph) + exp(-2*I*ph)*Ynm(1, 1, th, ph))/2) assert Znm(0, 0, th, ph).expand(func=True) == 1/(2*sqrt(pi)) assert Znm(1, -1, th, ph).expand(func=True) == (sqrt(3)*I*sin(th)*exp(I*ph)/(4*sqrt(pi)) - sqrt(3)*I*sin(th)*exp(-I*ph)/(4*sqrt(pi))) assert Znm(1, 0, th, ph).expand(func=True) == sqrt(3)*cos(th)/(2*sqrt(pi)) assert Znm(1, 1, th, ph).expand(func=True) == (-sqrt(3)*sin(th)*exp(I*ph)/(4*sqrt(pi)) - sqrt(3)*sin(th)*exp(-I*ph)/(4*sqrt(pi))) assert Znm(2, -1, th, ph).expand(func=True) == (sqrt(15)*I*sin(th)*exp(I*ph)*cos(th)/(4*sqrt(pi)) - sqrt(15)*I*sin(th)*exp(-I*ph)*cos(th)/(4*sqrt(pi))) assert Znm(2, 0, th, ph).expand(func=True) == 3*sqrt(5)*cos(th)**2/(4*sqrt(pi)) - sqrt(5)/(4*sqrt(pi)) assert Znm(2, 1, th, ph).expand(func=True) == (-sqrt(15)*sin(th)*exp(I*ph)*cos(th)/(4*sqrt(pi)) - sqrt(15)*sin(th)*exp(-I*ph)*cos(th)/(4*sqrt(pi)))
def as_real_imag(self, deep=True, **hints): """ Returns this function as a complex coordinate. """ from sympy import cos, sin if self.args[0].is_real: if deep: hints['complex'] = False return (self.expand(deep, **hints), S.Zero) else: return (self, S.Zero) if deep: re, im = self.args[0].expand(deep, **hints).as_real_imag() else: re, im = self.args[0].as_real_imag() return (sinh(re)*cos(im), cosh(re)*sin(im))
def test_issue_4892a(): A, z = symbols('A z') c = Symbol('c', nonzero=True) P1 = -A*exp(-z) P2 = -A/(c*t)*(sin(x)**2 + cos(y)**2) h1 = -sin(x)**2 - cos(y)**2 h2 = -sin(x)**2 + sin(y)**2 - 1 # there is still some non-deterministic behavior in integrate # or trigsimp which permits one of the following assert integrate(c*(P2 - P1), t) in [ c*(-A*(-h1)*log(c*t)/c + A*t*exp(-z)), c*(-A*(-h2)*log(c*t)/c + A*t*exp(-z)), c*( A* h1 *log(c*t)/c + A*t*exp(-z)), c*( A* h2 *log(c*t)/c + A*t*exp(-z)), (A*c*t - A*(-h1)*log(t)*exp(z))*exp(-z), (A*c*t - A*(-h2)*log(t)*exp(z))*exp(-z), ]
def __init__(self, n, index): self.dim = n if index == 2: self.degree = 2 r = sqrt(3) / 6 data = [ (1.0, numpy.array([numpy.full(n, 2*r)])), (+r, _s(n, -1, r)), (-r, _s(n, +1, r)), ] else: assert index == 3 self.degree = 3 n2 = n // 2 if n % 2 == 0 else (n-1)//2 i_range = range(1, 2*n+1) pts = [[ [sqrt(fr(2, 3)) * cos((2*k-1)*i*pi / n) for i in i_range], [sqrt(fr(2, 3)) * sin((2*k-1)*i*pi / n) for i in i_range], ] for k in range(1, n2+1)] if n % 2 == 1: sqrt3pm = numpy.full(2*n, 1/sqrt(3)) sqrt3pm[1::2] *= -1 pts.append(sqrt3pm) pts = numpy.vstack(pts).T data = [(fr(1, 2*n), pts)] self.points, self.weights = untangle(data) reference_volume = 2**n self.weights *= reference_volume return
def __init__(self, n): self.weights = numpy.full(n, 2*sympy.pi/n) self.points = numpy.column_stack([ [sympy.cos(sympy.pi * sympy.Rational(2*k, n)) for k in range(n)], [sympy.sin(sympy.pi * sympy.Rational(2*k, n)) for k in range(n)], ]) self.degree = n - 1 return
def _gen4_1(): pts = 2 * sqrt(5) * numpy.array([ [cos(2*i*pi/5) for i in range(5)], [sin(2*i*pi/5) for i in range(5)], ]).T data = [ (fr(7, 10), numpy.array([[0, 0]])), (fr(3, 50), pts), ] return 4, data # The boolean tells whether the factor 2*pi is already in the weights
