Extras¶
Solutions¶
Arithmetic operators¶
Solution 1¶
>>> from sympy import Symbol, Add, Poly
>>> x = Symbol('x')
>>> Add(*[ x**i for i in xrange(0, 10+1) ])
10 9 8 7 6 5 4 3 2
x + x + x + x + x + x + x + x + x + x + 1
>>> sum([ x**i for i in xrange(0, 10+1) ])
10 9 8 7 6 5 4 3 2
x + x + x + x + x + x + x + x + x + x + 1
>>> Poly([1]*11, x).as_expr()
10 9 8 7 6 5 4 3 2
x + x + x + x + x + x + x + x + x + x + 1
>>> def build_poly_1(n):
... return Add(*[ x**i for i in xrange(0, n+1) ])
...
>>> build_poly_1(1)
x + 1
>>> build_poly_1(2)
2
x + x + 1
>>> build_poly_1(3)
3 2
x + x + x + 1
>>> def build_poly_2(n):
... if n > 0:
... return Add(*[ x**i for i in xrange(0, n+1) ])
... else:
... return Add(*[ x**i for i in xrange(0, n-1, -1) ])
...
>>> build_poly_2(1)
x + 1
>>> build_poly_2(0)
1
>>> build_poly_2(-1)
1
1 + ─
x
Solution 2¶
>>> from sympy import Symbol
>>> x = Symbol('x')
>>> def nested_power(expr, n):
... if not n:
... return expr
... else:
... return expr**nested_power(expr, n-1)
...
>>> nested_power(x, 1)
x
x
>>> nested_power(x, 2)
⎛ x⎞
⎝x ⎠
x
>>> nested_power(x, 3)
⎛ ⎛ x⎞⎞
⎜ ⎝x ⎠⎟
⎝x ⎠
x
Building blocks of expressions¶
Solution 1¶
>>> from sympy import var, exp, sin, Integral, Eq
>>> var('a,b,c,x,C')
(a, b, c, x, C)
>>> expressions = [x**n, a*x**2 + b*x + c, 1/x, 1/(x**2 + 1), exp(a*x), sin(a*x) + b]
>>> expressions
⎡ n 2 1 1 a⋅x ⎤
⎢x , a⋅x + b⋅x + c, ─, ──────, ℯ , b + sin(a⋅x)⎥
⎢ x 2 ⎥
⎣ x + 1 ⎦
>>> integrals_table = []
>>> for expr in expressions:
... integral = Integral(expr, x)
... integrals_table.append((integral, integral.doit() + C))
...
>>> for integral, antiderivative in integrals_table:
... pprint(Eq(integral, antiderivative))
... print
...
⌠ n + 1
⎮ n x
⎮ x dx = C + ──────
⌡ n + 1
⌠ 3 2
⎮ ⎛ 2 ⎞ a⋅x b⋅x
⎮ ⎝a⋅x + b⋅x + c⎠ dx = C + ──── + ──── + c⋅x
⌡ 3 2
⌠
⎮ 1
⎮ ─ dx = C + log(x)
⎮ x
⌡
⌠
⎮ 1
⎮ ────── dx = C + atan(x)
⎮ 2
⎮ x + 1
⌡
⌠ a⋅x
⎮ a⋅x ℯ
⎮ ℯ dx = C + ────
⌡ a
⌠ cos(a⋅x)
⎮ (b + sin(a⋅x)) dx = C + b⋅x - ────────
⌡ a
⌠ n + 1
⎮ n x
⎮ x dx = C + ──────
⌡ n + 1
⌠ 3 2
⎮ ⎛ 2 ⎞ a⋅x b⋅x
⎮ ⎝a⋅x + b⋅x + c⎠ dx = C + ──── + ──── + c⋅x
⌡ 3 2
⌠
⎮ 1
⎮ ─ dx = C + log(x)
⎮ x
⌡
⌠
⎮ 1
⎮ ────── dx = C + atan(x)
⎮ 2
⎮ x + 1
⌡
⌠ a⋅x
⎮ a⋅x ℯ
⎮ ℯ dx = C + ────
⌡ a
⌠ cos(a⋅x)
⎮ (b + sin(a⋅x)) dx = C + b⋅x - ────────
⌡ a
Foreign types in SymPy¶
Solution 1¶
>>> from sympy.core.sympify import sympify, converter
>>> from gmpy import mpq
>>> def mpq_to_Rational(obj):
... return Rational(obj.numer(), obj.denom())
...
>>> converter[type(mpq(1))] = mpq_to_Rational
>>> sympify(mpq(1, 2))
1/2
>>> type(_)
<class 'sympy.core.numbers.Rational'>
Solution 2¶
>>> from sympy.core.sympify import converter, sympify, SympifyError
>>> from numpy import array, ndarray
>>> def ndarray_to_Tuple(obj):
... if len(obj.shape) == 1:
... return Tuple(*obj)
... else:
... raise SympifyError("only row NumPy arrays are allowed")
...
>>> converter[ndarray] = ndarray_to_Tuple
>>> sympify(array([1, 2, 3]))
Tuple(1, 2, 3)
>>> sympify(array([[1], [2], [3]]))
Traceback (most recent call last):
...
SympifyError: SympifyError: 'only row NumPy arrays are allowed'
The role of symbols¶
Solution 1¶
>>> from sympy import Add, Symbol, symbols, numbered_symbols
>>> def build_expression_1(name, n):
... return Add(*[ Symbol('%s_%d' % (name, i))**i for i in xrange(1, n+1) ])
...
>>> def build_expression_2(name, n):
... X = symbols('%s1:%d' % (name, n+1))
... return Add(*[ x**i for x, i in zip(X, xrange(1, n+1)) ])
...
>>> def build_expression_3(name, n):
... X = numbered_symbols(name, start=1)
... return Add(*[ x**i for x, i in zip(X, xrange(1, n+1)) ])
...
>>> build_expression_1('x', 5):
2 3 4 5
x₁ + x₂ + x₃ + x₄ + x₅
>>> build_expression_2('y', 5):
2 3 4 5
y₁ + y₂ + y₃ + y₄ + y₅
>>> build_expression_2('z', 5):
2 3 4 5
z₁ + z₂ + z₃ + z₄ + z₅
Obtaining parts of expressions¶
Solution 1¶
>>> from sympy.core import Atom, sympify
>>> def depth_1(expr):
... expr = sympify(expr)
...
... if isinstance(expr, Atom):
... return 1
... else:
... return 1 + max([ depth_1(arg) for arg in expr.args ])
...
>>> depth_1(117)
1
>>> def depth_2(expr):
... def _depth(expr):
... if isinstance(expr, Atom):
... return 1
... else:
... return 1 + max([ _depth(arg) for arg in expr.args ])
...
... return _depth(sympify(expr))
...
>>> depth_2(117)
1
Solution 2¶
>>> from sympy.core import Atom
>>> def depth(expr):
... if isinstance(expr, Atom):
... return 1
... else:
... if isinstance(expr, Basic):
... args = expr.args
... elif hasattr(expr, '__iter__'):
... args = expr
... else:
... raise ValueError("can't compute the depth of %s" % expr)
...
... return 1 + max([ depth(arg) for arg in args ])
...
>>> depth(x)
1
>>> depth(x + y)
2
>>> depth([x, y])
2
>>> depth((x, sin(x)))
3
>>> depth(set([(x, sin(x), sin(cos(x)))]))
4
Immutability of expressions¶
Solution 1¶
>>> var('x')
x
>>> x**x
x
x
>>> _.subs(x, x**x)
⎛ x⎞
⎝x ⎠
⎛ x⎞
⎝x ⎠
>>> _.subs(x, x**x)
⎛ ⎛ x⎞⎞
⎜ ⎝x ⎠⎟
⎜⎛ x⎞ ⎟
⎝⎝x ⎠ ⎠
⎛ ⎛ x⎞⎞
⎜ ⎝x ⎠⎟
⎜⎛ x⎞ ⎟
⎝⎝x ⎠ ⎠
>>> _.subs(x, x**x)
⎛ ⎛ ⎛ x⎞⎞⎞
⎜ ⎜ ⎝x ⎠⎟⎟
⎜ ⎜⎛ x⎞ ⎟⎟
⎜ ⎝⎝x ⎠ ⎠⎟
⎜⎛ ⎛ x⎞⎞ ⎟
⎜⎜ ⎝x ⎠⎟ ⎟
⎜⎜⎛ x⎞ ⎟ ⎟
⎝⎝⎝x ⎠ ⎠ ⎠
⎛ ⎛ ⎛ x⎞⎞⎞
⎜ ⎜ ⎝x ⎠⎟⎟
⎜ ⎜⎛ x⎞ ⎟⎟
⎜ ⎝⎝x ⎠ ⎠⎟
⎜⎛ ⎛ x⎞⎞ ⎟
⎜⎜ ⎝x ⎠⎟ ⎟
⎜⎜⎛ x⎞ ⎟ ⎟
⎝⎝⎝x ⎠ ⎠ ⎠
>>> expr = x**x
>>> for i in xrange(3):
... expr = (expr.base**expr.base)**(expr.exp**expr.exp)
...
>>> expr
⎛ ⎛ ⎛ x⎞⎞⎞
⎜ ⎜ ⎝x ⎠⎟⎟
⎜ ⎜⎛ x⎞ ⎟⎟
⎜ ⎝⎝x ⎠ ⎠⎟
⎜⎛ ⎛ x⎞⎞ ⎟
⎜⎜ ⎝x ⎠⎟ ⎟
⎜⎜⎛ x⎞ ⎟ ⎟
⎝⎝⎝x ⎠ ⎠ ⎠
⎛ ⎛ ⎛ x⎞⎞⎞
⎜ ⎜ ⎝x ⎠⎟⎟
⎜ ⎜⎛ x⎞ ⎟⎟
⎜ ⎝⎝x ⎠ ⎠⎟
⎜⎛ ⎛ x⎞⎞ ⎟
⎜⎜ ⎝x ⎠⎟ ⎟
⎜⎜⎛ x⎞ ⎟ ⎟
⎝⎝⎝x ⎠ ⎠ ⎠
Customizing built-in printers¶
Solution 1¶
>>> from sympy import Symbol, Lambda, sin
>>> from sympy.printing.pretty.pretty import PrettyPrinter
>>> from sympy.printing.pretty.stringpict import prettyForm
>>> x = Symbol('x')
>>> class LambdaPrettyPrinter(PrettyPrinter):
... def _print_Lambda(self, expr):
... args = expr.args[0].args + (expr.args[1],)
... pform_head = prettyForm('Lambda')
... pform_tail = self._print_seq(args, '(', ')')
... pform = prettyForm(*pform_head.right(pform_tail))
... return pform
...
>>> LambdaPrettyPrinter().doprint(Lambda(x, sin(1/x)))
⎛ ⎛1⎞⎞
Lambda⎜x, sin⎜─⎟⎟
⎝ ⎝x⎠⎠
Solution 2¶
>>> from sympy.printing.pretty.pretty import PrettyPrinter
>>> from sympy.printing.pretty.stringpict import prettyForm, stringPict
>>> class PolyPrettyPrinter(PrettyPrinter):
... def _print_Poly(self, poly):
... expr = poly.as_expr()
... gens = list(poly.gens)
... domain = poly.get_domain()
...
... pform = self._print(expr)
... pform = prettyForm(*stringPict.next(pform, ", "))
...
... for gen in gens:
... pform = prettyForm(*stringPict.next(pform, self._print(gen)))
... pform = prettyForm(*stringPict.next(pform, ", "))
...
... pform = prettyForm(*stringPict.next(pform, "domain="))
... pform = prettyForm(*stringPict.next(pform, self._print(domain)))
...
... pform = prettyForm(*pform.parens("(", ")", ifascii_nougly=True))
... pform = prettyForm(*prettyForm('Poly').right(pform))
...
... return pform
...
>>> PolyPrettyPrinter().doprint(Poly(x**2 + 1))
⎛ 2 ⎞
Poly⎝x + 1, x, domain=ℤ⎠
>>> PolyPrettyPrinter().doprint(Poly(x**2 + y/2))
⎛ 2 y ⎞
Poly⎜x + ─, x, y, domain=ℚ⎟
⎝ 2 ⎠
Implementing printers from scratch¶
Solution 1¶
>>> class ExtendedMathematicaPrinter(MathematicaPrinter):
... def _print_Pi(self, expr):
... return 'Pi'
...
>>> MathematicaPrinter().doprint(pi)
pi
>>> ExtendedMathematicaPrinter().doprint(pi)
Pi
Solution 2¶
>>> class ExtendedMathematicaPrinter(MathematicaPrinter):
... def _print_Mul(self, expr):
... prec = precedence(expr)
... return "*".join([ self.parenthesize(arg, prec) for arg in expr.args ])
...
>>> MathematicaPrinter().doprint(2*sin(x))
2*sin(x)
>>> ExtendedMathematicaPrinter().doprint(2*sin(x))
2*Sin[x]
Code generation¶
Solution 1¶
>>> from sympy.utilities.codegen import codegen
>>> from sympy import Symbol, chebyshevt
>>> x = Symbol('x')
>>> print codegen(("chebyshevt_20", chebyshevt(20, x)), "C", "file")[0][1]
/******************************************************************************
* Code generated with sympy 0.7.1 *
* *
* See http://www.sympy.org/ for more information. *
* *
* This file is part of 'project' *
******************************************************************************/
#include "file.h"
#include <math.h>
double chebyshevt_20(double x) {
return 524288*pow(x, 20) - 2621440*pow(x, 18) + 5570560*pow(x, 16) - 6553600*pow(x, 14) + 4659200*pow(x, 12) - 2050048*pow(x, 10) + 549120*pow(x, 8) - 84480*pow(x, 6) + 6600*pow(x, 4) - 200*pow(x, 2) + 1;
}
Solution 2¶
>>> from sympy.utilities.codegen import codegen
>>> from sympy import Symbol, chebyshevt
>>> x = Symbol('x')
>>> functions = [ ("chebyshevt_%d" % i, chebyshevt(i, x)) for i in xrange(0, 10) ]
>>> print codegen(functions, "F95", "file")[0][1]
!******************************************************************************
!* Code generated with sympy 0.7.1 *
!* *
!* See http://www.sympy.org/ for more information. *
!* *
!* This file is part of 'project' *
!******************************************************************************
REAL*8 function chebyshevt_0()
implicit none
chebyshevt_0 = 1
end function
REAL*8 function chebyshevt_1(x)
implicit none
REAL*8, intent(in) :: x
chebyshevt_1 = x
end function
REAL*8 function chebyshevt_2(x)
implicit none
REAL*8, intent(in) :: x
chebyshevt_2 = 2*x**2 - 1
end function
REAL*8 function chebyshevt_3(x)
implicit none
REAL*8, intent(in) :: x
chebyshevt_3 = 4*x**3 - 3*x
end function
REAL*8 function chebyshevt_4(x)
implicit none
REAL*8, intent(in) :: x
chebyshevt_4 = 8*x**4 - 8*x**2 + 1
end function
REAL*8 function chebyshevt_5(x)
implicit none
REAL*8, intent(in) :: x
chebyshevt_5 = 16*x**5 - 20*x**3 + 5*x
end function
REAL*8 function chebyshevt_6(x)
implicit none
REAL*8, intent(in) :: x
chebyshevt_6 = 32*x**6 - 48*x**4 + 18*x**2 - 1
end function
REAL*8 function chebyshevt_7(x)
implicit none
REAL*8, intent(in) :: x
chebyshevt_7 = 64*x**7 - 112*x**5 + 56*x**3 - 7*x
end function
REAL*8 function chebyshevt_8(x)
implicit none
REAL*8, intent(in) :: x
chebyshevt_8 = 128*x**8 - 256*x**6 + 160*x**4 - 32*x**2 + 1
end function
REAL*8 function chebyshevt_9(x)
implicit none
REAL*8, intent(in) :: x
chebyshevt_9 = 256*x**9 - 576*x**7 + 432*x**5 - 120*x**3 + 9*x
end function
Partial fraction decomposition¶
Solution 1¶
\(\frac{3 x + 5}{(2 x + 1)^2}\):
>>> f = (3*x + 5)/(2*x + 1)**2;f 3⋅x + 5 ────────── 2 (2⋅x + 1) >>> var('A:B') (A, B) >>> p1 = A/(2*x + 1) >>> p2 = B/(2*x + 1)**2 >>> p1, p2 ⎛ A B ⎞ ⎜───────, ──────────⎟ ⎜2⋅x + 1 2⎟ ⎝ (2⋅x + 1) ⎠ >>> h = sum((p1, p2));h A B ─────── + ────────── 2⋅x + 1 2 (2⋅x + 1) >>> hcombined = factor(together(h));hcombined 2⋅A⋅x + A + B ───────────── 2 (2⋅x + 1) >>> eq = Eq(hcombined, f);eq 2⋅A⋅x + A + B 3⋅x + 5 ───────────── = ────────── 2 2 (2⋅x + 1) (2⋅x + 1) >>> coeffeq = Eq(numer(eq.lhs), numer(eq.rhs));eq 2⋅A⋅x + A + B = 3⋅x + 5 >>> sol = solve_undetermined_coeffs(coeffeq, [A, B], x);sol {A: 3/2, B: 7/2} >>> solution = h.subs(sol);solution 3 7 ─────────── + ──────────── 2⋅(2⋅x + 1) 2 2⋅(2⋅x + 1) >>> apart(f, x) 3 7 ─────────── + ──────────── 2⋅(2⋅x + 1) 2 2⋅(2⋅x + 1) >>> solution == apart(f, x) True
\(\frac{3 x + 5}{(u x + v)^2}\):
>>> var('u v') (u, v) >>> f = (3*x + 5)/(u*x + v)**2;f 3⋅x + 5 ────────── 2 (u⋅x + v) >>> var('A:B') (A, B) >>> p1 = A/(u*x + v) >>> p2 = B/(u*x + v)**2 >>> p1, p2 ⎛ A B ⎞ ⎜───────, ──────────⎟ ⎜u⋅x + v 2⎟ ⎝ (u⋅x + v) ⎠ >>> h = sum((p1, p2));h A B ─────── + ────────── u⋅x + v 2 (u⋅x + v) >>> hcombined = factor(together(h));hcombined A⋅u⋅x + A⋅v + B ─────────────── 2 (u⋅x + v) >>> eq = Eq(hcombined, f);eq A⋅u⋅x + A⋅v + B 3⋅x + 5 ─────────────── = ────────── 2 2 (u⋅x + v) (u⋅x + v) >>> coeffeq = Eq(numer(eq.lhs), numer(eq.rhs));eq A⋅u⋅x + A⋅v + B 3⋅x + 5 ─────────────── = ────────── 2 2 (u⋅x + v) (u⋅x + v) >>> lhs, rhs = Poly(coeffeq.lhs, x), Poly(coeffeq.rhs, x) >>> lhs Poly(A*u*x + A*v + B, x, domain='ZZ[u,v,A,B]') >>> rhs Poly(3*x + 5, x, domain='ZZ') >>> coeffeqs = [ Eq(lhs.nth(i), rhs.nth(i)) for i in xrange(2) ];coeffeqs [A⋅v + B = 5, A⋅u = 3] >>> solution = h.subs(sol);solution 5⋅u - 3⋅v 3 ──────────── + ─────────── 2 u⋅(u⋅x + v) u⋅(u⋅x + v) >>> apart(f, x) # Note, you must include x! 5⋅u - 3⋅v 3 ──────────── + ─────────── 2 u⋅(u⋅x + v) u⋅(u⋅x + v) >>> solution == apart(f, x) True
\(\frac{(3 x + 5)^2}{(2 x + 1)^2}\):
>>> f = (3*x + 5)**2/(2*x + 1)**2;f 2 (3⋅x + 5) ────────── 2 (2⋅x + 1) >>> # Here, deg(p) == deg(q), so we have to use div() >>> p, q = numer(f), denom(f) >>> d, r = div(p, q) # p = d*q + r, i.e., f = p/q = d + r/q >>> d, r (9/4, 21⋅x + 91/4) >>> solution = d # This is the polynomial part of the solution >>> var('A:B') (A, B) >>> p1 = A/(2*x + 1) >>> p2 = B/(2*x + 1)**2 >>> p1, p2 ⎛ A B ⎞ ⎜───────, ──────────⎟ ⎜2⋅x + 1 2⎟ ⎝ (2⋅x + 1) ⎠ >>> h = sum((p1, p2));h A B ─────── + ────────── 2⋅x + 1 2 (2⋅x + 1) >>> hcombined = factor(together(h));hcombined 2⋅A⋅x + A + B ───────────── 2 (2⋅x + 1) >>> eq = Eq(hcombined, r/q);eq 2⋅A⋅x + A + B 21⋅x + 91/4 ───────────── = ─────────── 2 2 (2⋅x + 1) (2⋅x + 1) >>> coeffeq = Eq(numer(eq.lhs), r);eq # No need to call numer() here; we know it's r 2⋅A⋅x + A + B 21⋅x + 91/4 ───────────── = ─────────── 2 2 (2⋅x + 1) (2⋅x + 1) >>> coeffeqs = collect(coeffeq.lhs - coeffeq.rhs, x, evaluate=False);sol {A: 21/2, B: 49/4} >>> # Note, since we're just subtracting coeffeq.lhs from coeffeq.rhs, we could have >>> # done it without using Eq() in the first place. >>> sol = solve(coeffeqs.values()) >>> solution += h.subs(sol);solution # Note the += 9 21 49 ─ + ─────────── + ──────────── 4 2⋅(2⋅x + 1) 2 4⋅(2⋅x + 1) >>> apart(f, x) 9 21 49 ─ + ─────────── + ──────────── 4 2⋅(2⋅x + 1) 2 4⋅(2⋅x + 1) >>> solution == apart(f, x) True
Solution 2¶
Even though you can use Expr.coeff() to get the coefficient of \(x^n\) in a polynomial for \(n > 0\), it will not work for \(n = 0\). This is because Expr.coeff() considers any term with no numerical coefficient to be the coefficient of 1. From the docstring:
You can select terms with no rational coefficient:
>>> (x+2*y).coeff(1)
x
>>> (3+2*x+4*x**2).coeff(1)
Deriving trigonometric identities¶
Solution 1¶
>>> var('a,b')
(a, b)
>>> f = sin(a + b);f
sin(a + b)
>>> g = f.series(a, 0, 10)
>>> g
2 3 4 5 6 7 8 9
a ⋅sin(b) a ⋅cos(b) a ⋅sin(b) a ⋅cos(b) a ⋅sin(b) a ⋅cos(b) a ⋅sin(b) a ⋅cos(b)
sin(b) + a⋅cos(b) - ───────── - ───────── + ───────── + ───────── - ───────── - ───────── + ───────── + ───────── + O(a**10)
2 6 24 120 720 5040 40320 362880
>>> g = collect(g, [sin(b), cos(b)]).removeO()
>>> g
⎛ 8 6 4 2 ⎞ ⎛ 9 7 5 3 ⎞
⎜ a a a a ⎟ ⎜ a a a a ⎟
⎜───── - ─── + ── - ── + 1⎟⋅sin(b) + ⎜────── - ──── + ─── - ── + a⎟⋅cos(b)
⎝40320 720 24 2 ⎠ ⎝362880 5040 120 6 ⎠
>>> sina = sin(a).series(a, 0, 10).removeO()
>>> sina
9 7 5 3
a a a a
────── - ──── + ─── - ── + a
362880 5040 120 6
>>> h = g.subs(sina, sin(a))
>>> h
⎛ 8 6 4 2 ⎞
⎜ a a a a ⎟
⎜───── - ─── + ── - ── + 1⎟⋅sin(b) + sin(a)⋅cos(b)
⎝40320 720 24 2 ⎠
>>> cosa = cos(a).series(a, 0, 10).removeO()
>>> cosa
8 6 4 2
a a a a
───── - ─── + ── - ── + 1
40320 720 24 2
>>> h = h.subs(cosa, cos(a))
>>> h
sin(a)⋅cos(b) + sin(b)⋅cos(a)
>>> eq = Eq(f, h)
>>> eq
sin(a + b) = sin(a)⋅cos(b) + sin(b)⋅cos(a)
>>> eq.lhs.expand(trig=True) == eq.rhs
True
Solution 2¶
>>> var('a,b')
(a, b)
>>> f = cos(a + b);f
cos(a + b)
>>> g = f.series(a, 0, 10)
>>> g
2 3 4 5 6 7 8 9
a ⋅cos(b) a ⋅sin(b) a ⋅cos(b) a ⋅sin(b) a ⋅cos(b) a ⋅sin(b) a ⋅cos(b) a ⋅sin(b)
cos(b) - a⋅sin(b) - ───────── + ───────── + ───────── - ───────── - ───────── + ───────── + ───────── - ───────── + O(a**10)
2 6 24 120 720 5040 40320 362880
>>> g = collect(g, [sin(b), cos(b)]).removeO()
>>> g
⎛ 8 6 4 2 ⎞ ⎛ 9 7 5 3 ⎞
⎜ a a a a ⎟ ⎜ a a a a ⎟
⎜───── - ─── + ── - ── + 1⎟⋅cos(b) + ⎜- ────── + ──── - ─── + ── - a⎟⋅sin(b)
⎝40320 720 24 2 ⎠ ⎝ 362880 5040 120 6 ⎠
>>> sina = sin(a).series(a, 0, 10).removeO()
>>> sina
9 7 5 3
a a a a
────── - ──── + ─── - ── + a
362880 5040 120 6
>>> # Note that subs will not work directly with sina, because -sina appears
>>> # in g. We get around it by substituting -sina.
>>> h = g.subs(-sina, -sin(a))
>>> h
⎛ 8 6 4 2 ⎞
⎜ a a a a ⎟
⎜───── - ─── + ── - ── + 1⎟⋅cos(b) - sin(a)⋅sin(b)
⎝40320 720 24 2 ⎠
>>> cosa = cos(a).series(a, 0, 10).removeO()
>>> cosa
8 6 4 2
a a a a
───── - ─── + ── - ── + 1
40320 720 24 2
>>> h = h.subs(cosa, cos(a))
>>> h
-sin(a)⋅sin(b) + cos(a)⋅cos(b)
>>> eq = Eq(f, h)
>>> eq
cos(a + b) = -sin(a)⋅sin(b) + cos(a)⋅cos(b)
>>> eq.lhs.expand(trig=True) == eq.rhs
True
Not only symbolics: numerical computing¶
Solution 1¶
>>> f = x**(1 - log(log(log(log(1/x)))))
>>> f.subs(x, pi).evalf()
0.7302900194855046292629153268423316194247957512439873318 - 1.308036181902133609036447992994107130426900293406430474⋅ⅈ
Solution 2¶
>>> E_million_digits = str(E.evalf(n=1000000))
>>> pi_million_digits = str(pi.evalf(n=1000000))
>>> # We can use the .find() method of str.
>>> # First, let's see what we should offset the value
>>> # (it's not immediately obvious because of 0 indexing
>>> # and because of the '.' in the string.
>>> E_million_digits[:5]
2.718
>>> E_million_digits.find('71')
2
>>> # So, counting the 2 as the first digit, we see that
>>> # there is no need for any offset!
>>> E999999 = E_million_digits.find('999999')
>>> E999999
384341
>>> "...%s..." % E_million_digits[E999999:E999999+10]
...9999999951...
>>> # We see that there are actually 8 consecutive 9's here
>>> pi789 = pi_million_digits.find('789')
>>> pi789
352
>>> "...%s..." % pi_million_digits[pi789:pi789+10]
...7892590360...
>>> # Now, double check the values from Dinosaur Comics.
>>> E789 = E_million_digits.find('789')
>>> E789
2501
>>> "...%s..." % E_million_digits[E789:E789+10]
...7893652605...
>>> pi999999 = pi_million_digits.find('999999')
>>> pi999999
763
>>> # Hey, this one was wrong in the webcomic (T-Rex said "762nd")!
>>> "...%s..." % pi_million_digits[pi999999:pi999999+10]
...9999998372...
Solution 3¶
>>> limit((erf(x - exp(-exp(x))) - erf(x))*exp(exp(x))*exp(x**2), x, oo)
-2
─────
⎽⎽⎽
╲╱ π
Summing roots of polynomials¶
Solution 1¶
\(f = z^5 + z + a\) and \(g = \frac{1}{z + 1}\):
>>> var('a') a >>> f = z**5 + z + a >>> f 5 a + z + z >>> g = 1/(z + 1) >>> g 1 ───── z + 1 >>> R = var('r:5') >>> G = together(sum( [ g.subs(z, r) for r in R] )) >>> G (r₀ + 1)⋅(r₁ + 1)⋅(r₂ + 1)⋅(r₃ + 1) + (r₀ + 1)⋅(r₁ + 1)⋅(r₂ + 1)⋅(r₄ + 1) + (r₀ + 1)⋅(r₁ + 1)⋅(r₃ + 1)⋅(r₄ + 1) + (r₀ + 1)⋅(r₂ + 1)⋅(r₃ + 1)⋅(r₄ + 1) + (r₁ + 1)⋅(r₂ + 1)⋅(r₃ + 1)⋅(r₄ + 1) ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── (r₀ + 1)⋅(r₁ + 1)⋅(r₂ + 1)⋅(r₃ + 1)⋅(r₄ + 1) >>> V = viete(f, R, z) >>> for lhs, rhs in V: ... pprint(Eq(lhs, rhs)) ... r₀ + r₁ + r₂ + r₃ + r₄ = 0 r₀⋅r₁ + r₀⋅r₂ + r₀⋅r₃ + r₀⋅r₄ + r₁⋅r₂ + r₁⋅r₃ + r₁⋅r₄ + r₂⋅r₃ + r₂⋅r₄ + r₃⋅r₄ = 0 r₀⋅r₁⋅r₂ + r₀⋅r₁⋅r₃ + r₀⋅r₁⋅r₄ + r₀⋅r₂⋅r₃ + r₀⋅r₂⋅r₄ + r₀⋅r₃⋅r₄ + r₁⋅r₂⋅r₃ + r₁⋅r₂⋅r₄ + r₁⋅r₃⋅r₄ + r₂⋅r₃⋅r₄ = 0 r₀⋅r₁⋅r₂⋅r₃ + r₀⋅r₁⋅r₂⋅r₄ + r₀⋅r₁⋅r₃⋅r₄ + r₀⋅r₂⋅r₃⋅r₄ + r₁⋅r₂⋅r₃⋅r₄ = 1 r₀⋅r₁⋅r₂⋅r₃⋅r₄ = -a >>> p = expand(numer(G)) >>> q = expand(denom(G)) >>> (P, Q), mapping = symmetrize((p, q), R, formal=True) >>> P (4⋅s₁ + 3⋅s₂ + 2⋅s₃ + s₄ + 5, 0) >>> Q (s₁ + s₂ + s₃ + s₄ + s₅ + 1, 0) >>> for s, poly in mapping: ... pprint(Eq(s, poly)) ... s₁ = r₀ + r₁ + r₂ + r₃ + r₄ s₂ = r₀⋅r₁ + r₀⋅r₂ + r₀⋅r₃ + r₀⋅r₄ + r₁⋅r₂ + r₁⋅r₃ + r₁⋅r₄ + r₂⋅r₃ + r₂⋅r₄ + r₃⋅r₄ s₃ = r₀⋅r₁⋅r₂ + r₀⋅r₁⋅r₃ + r₀⋅r₁⋅r₄ + r₀⋅r₂⋅r₃ + r₀⋅r₂⋅r₄ + r₀⋅r₃⋅r₄ + r₁⋅r₂⋅r₃ + r₁⋅r₂⋅r₄ + r₁⋅r₃⋅r₄ + r₂⋅r₃⋅r₄ s₄ = r₀⋅r₁⋅r₂⋅r₃ + r₀⋅r₁⋅r₂⋅r₄ + r₀⋅r₁⋅r₃⋅r₄ + r₀⋅r₂⋅r₃⋅r₄ + r₁⋅r₂⋅r₃⋅r₄ s₅ = r₀⋅r₁⋅r₂⋅r₃⋅r₄ >>> subslist = [ (s, c) for (s, _), (_, c) in zip(mapping, V) ] >>> subslist [(s₁, 0), (s₂, 0), (s₃, 0), (s₄, 1), (s₅, -a)] >>> P[0].subs(subslist)/Q[0].subs(subslist) 6 ────── -a + 2 >>> # Note, we must give a variable to RootSum, as f has more than one Symbol >>> cancel(P[0].subs(subslist)/Q[0].subs(subslist) - RootSum(f, Lambda(z, g), z)) == 0 True
\(f = z^5 + z + a\) and \(g = \frac{1}{z + b}\):
>>> var('a b') (a, b) >>> f = z**5 + z + a >>> f 5 a + z + z >>> g = 1/(z + b) >>> g 1 ───── b + z >>> R = var('r:5') >>> G = together(sum( [ g.subs(z, r) for r in R] )) >>> G (b + r₀)⋅(b + r₁)⋅(b + r₂)⋅(b + r₃) + (b + r₀)⋅(b + r₁)⋅(b + r₂)⋅(b + r₄) + (b + r₀)⋅(b + r₁)⋅(b + r₃)⋅(b + r₄) + (b + r₀)⋅(b + r₂)⋅(b + r₃)⋅(b + r₄) + (b + r₁)⋅(b + r₂)⋅(b + r₃)⋅(b + r₄) ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── (b + r₀)⋅(b + r₁)⋅(b + r₂)⋅(b + r₃)⋅(b + r₄) >>> V = viete(f, R, z) >>> for lhs, rhs in V: ... pprint(Eq(lhs, rhs)) ... r₀ + r₁ + r₂ + r₃ + r₄ = 0 r₀⋅r₁ + r₀⋅r₂ + r₀⋅r₃ + r₀⋅r₄ + r₁⋅r₂ + r₁⋅r₃ + r₁⋅r₄ + r₂⋅r₃ + r₂⋅r₄ + r₃⋅r₄ = 0 r₀⋅r₁⋅r₂ + r₀⋅r₁⋅r₃ + r₀⋅r₁⋅r₄ + r₀⋅r₂⋅r₃ + r₀⋅r₂⋅r₄ + r₀⋅r₃⋅r₄ + r₁⋅r₂⋅r₃ + r₁⋅r₂⋅r₄ + r₁⋅r₃⋅r₄ + r₂⋅r₃⋅r₄ = 0 r₀⋅r₁⋅r₂⋅r₃ + r₀⋅r₁⋅r₂⋅r₄ + r₀⋅r₁⋅r₃⋅r₄ + r₀⋅r₂⋅r₃⋅r₄ + r₁⋅r₂⋅r₃⋅r₄ = 1 r₀⋅r₁⋅r₂⋅r₃⋅r₄ = -a >>> p = expand(numer(G)) >>> q = expand(denom(G)) >>> (P, Q), mapping = symmetrize((p, q), R, formal=True) >>> P ⎛ 4 3 2 ⎞ ⎝5⋅b + 4⋅b ⋅s₁ + 3⋅b ⋅s₂ + 2⋅b⋅s₃ + s₄, 0⎠ >>> Q ⎛ 5 4 3 2 ⎞ ⎝b + b ⋅s₁ + b ⋅s₂ + b ⋅s₃ + b⋅s₄ + s₅, 0⎠ >>> for s, poly in mapping: ... pprint(Eq(s, poly)) ... s₁ = r₀ + r₁ + r₂ + r₃ + r₄ s₂ = r₀⋅r₁ + r₀⋅r₂ + r₀⋅r₃ + r₀⋅r₄ + r₁⋅r₂ + r₁⋅r₃ + r₁⋅r₄ + r₂⋅r₃ + r₂⋅r₄ + r₃⋅r₄ s₃ = r₀⋅r₁⋅r₂ + r₀⋅r₁⋅r₃ + r₀⋅r₁⋅r₄ + r₀⋅r₂⋅r₃ + r₀⋅r₂⋅r₄ + r₀⋅r₃⋅r₄ + r₁⋅r₂⋅r₃ + r₁⋅r₂⋅r₄ + r₁⋅r₃⋅r₄ + r₂⋅r₃⋅r₄ s₄ = r₀⋅r₁⋅r₂⋅r₃ + r₀⋅r₁⋅r₂⋅r₄ + r₀⋅r₁⋅r₃⋅r₄ + r₀⋅r₂⋅r₃⋅r₄ + r₁⋅r₂⋅r₃⋅r₄ s₅ = r₀⋅r₁⋅r₂⋅r₃⋅r₄ >>> subslist = [ (s, c) for (s, _), (_, c) in zip(mapping, V) ] >>> subslist [(s₁, 0), (s₂, 0), (s₃, 0), (s₄, 1), (s₅, -a)] >>> P[0].subs(subslist)/Q[0].subs(subslist) 4 5⋅b + 1 ─────────── 5 -a + b + b >>> # Note, we must give a variable to RootSum, as f has more than one Symbol >>> cancel(P[0].subs(subslist)/Q[0].subs(subslist) - RootSum(f, Lambda(z, g), z)) == 0 True
Solution 2¶
Table Of Contents
- Extras
- Solutions
- Arithmetic operators
- Building blocks of expressions
- Foreign types in SymPy
- The role of symbols
- Obtaining parts of expressions
- Immutability of expressions
- Turning strings into expressions
- Customizing built-in printers
- Implementing printers from scratch
- Code generation
- Partial fraction decomposition
- Deriving trigonometric identities
- Not only symbolics: numerical computing
- Summing roots of polynomials
- Vertex \(k\)–coloring of graphs
- Algebraic geometry
- Solutions