On 15 June 2013 10:53, Amit Saha <amitsaha...@gmail.com> wrote: > Would any of you have any teaching (or substantial self learning) > experience with a library for Symbolic math?
I have taught using Maple. I haven't with Sympy but I do use it myself. > I am currently exploring sympy (http://sympy.org) as part of writing a > book chapter and would like to know if there any better/easier option > out there which can successfully introduce symbolic math to young > programmers. In terms of free software sympy/sage are probably the best options. I guess that my preference for teaching symbolic computation would be to begin with isympy (sympy mode in the IPython console). I think that this is initially better than the notebook modes of sympy, Maple, etc. because it is clear what is input and what is output and what functions are being called. That way you can gain a sense of the nuts and bolts for how symbolic computation works without getting confused by the "magic" of most of these systems. Later most people would probably find it more practical to use the isympy notebooks or sage for serious work. If you use isympy with the -a flag then missing symbols are automatically created solving the inconvenience that Jim referred to. Also if you have a decent console you can pass '-p unicode' to get all the Greek letters (this is automatic on e.g. Ubuntu but not on Windows): $ isympy -p unicode -a IPython console for SymPy 0.7.2 (Python 2.7.5-32-bit) (ground types: python) These commands were executed: >>> from __future__ import division >>> from sympy import * >>> x, y, z, t = symbols('x y z t') >>> k, m, n = symbols('k m n', integer=True) >>> f, g, h = symbols('f g h', cls=Function) Documentation can be found at http://www.sympy.org In [1]: alpha + kappa Out[1]: α + κ Unfortunately the auto-symbols feature doesn't work so well when a particular symbol is unexpectedly defined e.g. the beta function: In [2]: alpha + beta --------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-2-92b8e9c3b3ff> in <module>() ----> 1 alpha + beta TypeError: unsupported operand type(s) for +: 'Symbol' and 'function' Otherwise I think that sympy has a good feature set. Some examples: Get the general formula for the sum of i**2 with i from 1 to n: In [1]: summation(i**2, (i, 1, n)) Out[1]: 3 2 n n n -- + -- + - 3 2 6 Get the sum of 2^(-i) with i from 0 to infinity: In [3]: summation(2**-i, (i, 0, oo)) Out[3]: 2 Quickly compute the derivative of some hairy function: In [5]: (x ** 2 * exp(x) * (sin(omega*x + phi))).diff(x) Out[5]: 2 x 2 x x omega*x *e *cos(omega*x + phi) + x *e *sin(omega*x + phi) + 2*x*e *sin(omega*x + phi) If you have matplotlib installed then it can quickly plot things for you. Otherwise it falls back to a great ascii poltting library: In [11]: plot(x**2) 100 | | . . | . . | . . | . . | . . | . . | . 50.0165 | --------.--------------------------------------..------ | . . | . . | . . | .. .. | . . | .. .. | .. .. | .. .. 0.03305 | .............. -10 0 10 Out[11]: Plot object containing: [0]: cartesian line: x**2 for x over (-10.0, 10.0) Sympy is good for simple interactive work. It is also suitable for complex problems but some of the interfaces are a little clunky when you start needing to manipulate things programmatically. Here's a non-trivial script that I used to compute the Butcher table for the 4th-order Gauss-Legendre method: #!/usr/bin/env python # # Compute the Butcher table for the Gauss-Legendre 4th order scheme from sympy import * c1, c2, a1, a2, b1, b2, k1, k2, xn, h = symbols('c1 c2 a1 a2 b1 b2 k1 k2 xn h') # Interpolation points # Roots of the Legendre polynomial (transformed) p1 = Rational(1, 2) - Rational(1, 6)*sqrt(3) p2 = Rational(1, 2) + Rational(1, 6)*sqrt(3) # Imagine that we are integrating the differential equation dxdt = x # In that case the solution to the implicit problem for the stages is given by # eq1 and eq2 below. eq1 = h * (xn + a1*k1 +a2*k2) - k1 eq2 = h * (xn + b1*k1 +b2*k2) - k2 # Lets solve the implicit problem to get k1 and k2 in terms of the table # coefficients. soln = solve([eq1, eq2], [k1, k2]) # True solution of the differential equation xnp1 = xn * exp(h) # Estimate that we get from our integration step xnp1bar = (xn + c1*soln[k1] + c2*soln[k2]) # The relative error between the true solution and estimate E = ((xnp1bar - xnp1) / xnp1) # Utility function to extract the terms from the Maclaurin series below # Must be an easier way to do this. def get_terms(expr, var): def prep(expr): expr = expr.subs(a2, p1 - a1).subs(b2, p2 - b1) expr = collect(cancel(expand(expr)), h) return expr results = {} order = -1 while order < 5: coeff, order = expr.leadterm(h) results[var**order] = prep(coeff) expr = cancel(expr - coeff * var ** order) return results # We'll expand the error as a Maclaurin series and try to cancel the first # three terms. Eh = get_terms(series(E, h, n=6), h) # Solve for the values of c, a1 and b1 that cancel the first three terms of # the series expansion. c_sol = solve([Eh[h**2], 1-c1-c2], [c1, c2]) for v in h**3, h**4: Eh[v] = cancel(expand(Eh[v].subs(c1, c_sol[c1]).subs(c2, c_sol[c2]))) (a_sol, b_sol), = solve([Eh[h**3], Eh[h**4]], [a1, b1]) # Construct as a Matrix for pretty-printing table = Matrix([ [p1 , a_sol, p1 - a_sol], [p2 , b_sol, p2 - b_sol], [nan, c_sol[c1], c_sol[c2]], ]) print('Butcher Table for Gauss-Legendre 4th order') pprint(table) The output (after ~5 seconds) is: $ python gl4.py Butcher Table for Gauss-Legendre 4th order [ ___ ___ ] [ \/ 3 1 \/ 3 1] [- ----- + - 1/4 - ----- + -] [ 6 2 6 4] [ ] [ ___ ___ ] [ \/ 3 1 1 \/ 3 ] [ ----- + - - + ----- 1/4 ] [ 6 2 4 6 ] [ ] [ nan 1/2 1/2 ] Compare this with the second table here: http://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_method Oscar _______________________________________________ Tutor maillist - Tutor@python.org To unsubscribe or change subscription options: http://mail.python.org/mailman/listinfo/tutor
