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 ⎠    ⎠          ⎠

Turning strings into expressions

Solution 1

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

Vertex \(k\)–coloring of graphs

Solution 1

Solution 2

Algebraic geometry

Solution 1