From: W. Trevor King <wking@tremily.us>
Date: Wed, 27 Feb 2013 12:51:53 +0000 (-0500)
Subject: posts:sympy: Add post demonstrating usage
X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=1cfc5217717767e745bf4a742c3ea5fd0dd0ca52;p=blog.git

posts:sympy: Add post demonstrating usage
---

diff --git a/posts/SymPy.mdwn_itex b/posts/SymPy.mdwn_itex
new file mode 100644
index 0000000..e8a7cde
--- /dev/null
+++ b/posts/SymPy.mdwn_itex
@@ -0,0 +1,121 @@
+[SymPy][] is a [[Python]] library for symbolic mathematics.  To give
+you a feel for how it works, lets extrapolate the extremum location
+for $f(x)$ given a quadratic model:
+
+\[
+  f(x) = A x^2 + B x + C
+\]
+
+and three known values:
+
+\[
+\begin{aligned}
+  f(a) &= A a^2 + B a + C  \\
+  f(b) &= A b^2 + B b + C  \\
+  f(c) &= A c^2 + B c + C
+\end{aligned}
+\]
+
+Rephrase as a matrix equation:
+
+\[
+\begin{pmatrix}
+  f(a)  \\
+  f(b)  \\
+  f(c)
+\end{pmatrix}
+  =
+\begin{pmatrix}
+  a^2  &  a  &  1  \\
+  b^2  &  b  &  1  \\
+  c^2  &  c  &  1
+\end{pmatrix}
+  \cdot
+\begin{pmatrix}
+  A  \\
+  B  \\
+  C
+\end{pmatrix}
+\]
+
+So the solutions for $A$, $B$, and $C$ are:
+
+\[
+\begin{pmatrix}
+  A  \\
+  B  \\
+  C
+\end{pmatrix}
+  =
+\begin{pmatrix}
+  a^2  &  a  &  1  \\
+  b^2  &  b  &  1  \\
+  c^2  &  c  &  1
+\end{pmatrix}^{-1}
+  \cdot
+\begin{pmatrix}
+  f(a)  \\
+  f(b)  \\
+  f(c)
+\end{pmatrix}
+  =
+\begin{pmatrix}
+  \text{long}  \\
+  \text{complicated}  \\
+  \text{stuff}
+\end{pmatrix}
+\]
+
+Now that we've found the model parameters, we need to find the $x$
+coordinate of the extremum.
+
+\[
+  \frac{\mathrm{d}f}{\mathrm{d}x} = 2 A x + B  \;,
+\]
+
+which is zero when
+
+\[
+\begin{aligned}
+  2 A x &= -B  \\
+  x &= \frac{-B}{2 A}
+\end{aligned}
+\]
+
+Here's the solution in SymPy:
+
+    >>> from sympy import Symbol, Matrix, factor, expand, pprint, preview
+    >>> a = Symbol('a')
+    >>> b = Symbol('b')
+    >>> c = Symbol('c')
+    >>> fa = Symbol('fa')
+    >>> fb = Symbol('fb')
+    >>> fc = Symbol('fc')
+    >>> M = Matrix([[a**2, a, 1], [b**2, b, 1], [c**2, c, 1]])
+    >>> F = Matrix([[fa],[fb],[fc]])
+    >>> ABC = M.inv() * F
+    >>> A = ABC[0,0]
+    >>> B = ABC[1,0]
+    >>> x = -B/(2*A)
+    >>> x = factor(expand(x))
+    >>> pprint(x)
+     2       2       2       2       2       2   
+    a *fb - a *fc - b *fa + b *fc + c *fa - c *fb
+    ---------------------------------------------
+     2*(a*fb - a*fc - b*fa + b*fc + c*fa - c*fb) 
+    >>> preview(x, viewer='pqiv')
+
+Where `pqiv` is the executable for [pqiv][], my preferred image
+viewer.  With a bit of additional factoring, that is:
+
+\[
+  x = \frac{a^2 [f(b)-f(c)] + b^2 [f(c) - f(a)] + c^2 [f(a) - f(b)]}
+           {2\cdot\{a [f(b) - f(c)] + b [f(c) - f(a)] + c [f(a) - f(b)]\}}
+\]
+
+[SymPy]: http://sympy.org/
+[pqiv]: http://www.pberndt.com/Programme/Linux/pqiv/
+
+[[!tag tags/python]]
+[[!tag tags/teaching]]
+[[!tag tags/tools]]