posts/SymPy.mdwn_itex

   1 [SymPy][] is a [[Python]] library for symbolic mathematics.  To give
   2 you a feel for how it works, lets extrapolate the extremum location
   3 for $f(x)$ given a quadratic model:
   4
   5 \[
   6   f(x) = A x^2 + B x + C
   7 \]
   8
   9 and three known values:
  10
  11 \[
  12 \begin{aligned}
  13   f(a) &= A a^2 + B a + C  \\
  14   f(b) &= A b^2 + B b + C  \\
  15   f(c) &= A c^2 + B c + C
  16 \end{aligned}
  17 \]
  18
  19 Rephrase as a matrix equation:
  20
  21 \[
  22 \begin{pmatrix}
  23   f(a)  \\
  24   f(b)  \\
  25   f(c)
  26 \end{pmatrix}
  27   =
  28 \begin{pmatrix}
  29   a^2  &  a  &  1  \\
  30   b^2  &  b  &  1  \\
  31   c^2  &  c  &  1
  32 \end{pmatrix}
  33   \cdot
  34 \begin{pmatrix}
  35   A  \\
  36   B  \\
  37   C
  38 \end{pmatrix}
  39 \]
  40
  41 So the solutions for $A$, $B$, and $C$ are:
  42
  43 \[
  44 \begin{pmatrix}
  45   A  \\
  46   B  \\
  47   C
  48 \end{pmatrix}
  49   =
  50 \begin{pmatrix}
  51   a^2  &  a  &  1  \\
  52   b^2  &  b  &  1  \\
  53   c^2  &  c  &  1
  54 \end{pmatrix}^{-1}
  55   \cdot
  56 \begin{pmatrix}
  57   f(a)  \\
  58   f(b)  \\
  59   f(c)
  60 \end{pmatrix}
  61   =
  62 \begin{pmatrix}
  63   \text{long}  \\
  64   \text{complicated}  \\
  65   \text{stuff}
  66 \end{pmatrix}
  67 \]
  68
  69 Now that we've found the model parameters, we need to find the $x$
  70 coordinate of the extremum.
  71
  72 \[
  73   \frac{\mathrm{d}f}{\mathrm{d}x} = 2 A x + B  \;,
  74 \]
  75
  76 which is zero when
  77
  78 \[
  79 \begin{aligned}
  80   2 A x &= -B  \\
  81   x &= \frac{-B}{2 A}
  82 \end{aligned}
  83 \]
  84
  85 Here's the solution in SymPy:
  86
  87     >>> from sympy import Symbol, Matrix, factor, expand, pprint, preview
  88     >>> a = Symbol('a')
  89     >>> b = Symbol('b')
  90     >>> c = Symbol('c')
  91     >>> fa = Symbol('fa')
  92     >>> fb = Symbol('fb')
  93     >>> fc = Symbol('fc')
  94     >>> M = Matrix([[a**2, a, 1], [b**2, b, 1], [c**2, c, 1]])
  95     >>> F = Matrix([[fa],[fb],[fc]])
  96     >>> ABC = M.inv() * F
  97     >>> A = ABC[0,0]
  98     >>> B = ABC[1,0]
  99     >>> x = -B/(2*A)
 100     >>> x = factor(expand(x))
 101     >>> pprint(x)
 102      2       2       2       2       2       2
 103     a *fb - a *fc - b *fa + b *fc + c *fa - c *fb
 104     ---------------------------------------------
 105      2*(a*fb - a*fc - b*fa + b*fc + c*fa - c*fb)
 106     >>> preview(x, viewer='pqiv')
 107
 108 Where `pqiv` is the executable for [pqiv][], my preferred image
 109 viewer.  With a bit of additional factoring, that is:
 110
 111 \[
 112   x = \frac{a^2 [f(b)-f(c)] + b^2 [f(c) - f(a)] + c^2 [f(a) - f(b)]}
 113            {2\cdot\{a [f(b) - f(c)] + b [f(c) - f(a)] + c [f(a) - f(b)]\}}
 114 \]
 115
 116 [SymPy]: http://sympy.org/
 117 [pqiv]: http://www.pberndt.com/Programme/Linux/pqiv/
 118
 119 [[!tag tags/python]]
 120 [[!tag tags/teaching]]
 121 [[!tag tags/tools]]