[Edu-sig] the plot question (in another sense)

[Kirby Urner]

At 01:28 PM 1/26/2003 -0800, Kirby Urner wrote:

>We could avoid the proofs for a bit and just play around with the first
>expansion, comparing outputs on both sides of the equals sign.  I.e. in
>Python we have math.exp(x) for raising e to the xth power.

Clearly if we're doing a polynomial expansion equal to the cosine of
x, its graph'll have to be the same as cosine's.  (x,cos(x)) crosses
the x-axis at all these roots (e.g. pi/2, 3pi/2, 5pi/2 etc.).  So the
polynomial needs another root, and therefore one higher degree, for
every new x-axis crossing.

The polynomial expansion for cosine is [1/0!, 0, -1/2!, 0, 1/4!, 0, -1/6!...]
where these are coefficients for x^0 x^1 x^2 etc. --except the terms with odd
exponents have zero coefficients and so go away.

Using the generator approach again:

  def costerms(x,t):
      n,result = 0,1
      maxterm = 2*t
      if n==0: yield 1
      while 1:
         n += 2
         if n >= maxterm: return   # bug fixed (>=) -- one too many terms 
before
         result *= -x*x/(n*(n-1))
         yield result

  def polycos(x,t):
      return reduce(add,[i for i in costerms(x,t)])

  >>> polycos(3,20)
  -0.98999249660044553
  >>> math.cos(3)
  -0.98999249660044542

So now I'd like see graphs of these ever higher degree polynomials that
criss-cross the x-axis, just like cosine does.  In J, I'd just go something 
like:

   load 'plot';'trig'
   polycos =: (cos t. i.20) & p.  NB. coefficients out 20 terms, including 
zeros
   domain  =: (i: 80) % 10        NB. from -8 to 8 in increments of 1/10
   plot domain; polycos domain    NB. do the graph

and I'd be looking at my wavy line on a graph in a window.

One approach I've used is to use a module I wrote for outputting to Povray,
the ray tracing engine.  I just output the (x,y) points as little spheres
and link them with cylinders.  This gives a graph with a pleasing 3d look
-- can be colorful too.

However, since I've recently added wxPython to Python 2.2 (and am using the
PyCrust that comes with it), tonight my experiment will by to use a package
I've not tried before: SciPy (VPython also has graphing capabilities BTW).
I've just grabbed the Windows binary from here:
http://www.scipy.org/site_content/download_list

It seems to install a *lot* of stuff, but then I see I still need Numeric,
which the VPython uninstaller recently stole off my machine (as Arthur
warned me it would).  So I'll go to http://www.pfdubois.com/numpy/ and grab
another copy... (Numpy is pretty small, only like 446K -- downloading *very*
slowly though, looks like I need to try again... ah, *much* faster from
Chapel Hill, NC).

OK, now scipy imports without errors.  Time to check the docs and see how
to do a quick plot... http://www.scipy.org/site_content/tutorials/plot_tutorial

Hmmmm... Strangely, when I used Vim to edit code, then reload my module
in PyCrust, I sometimes get error messages referring to old code that's
gone from the module.  I wonder if there's a namespace conflict with
scipy somehow.  This 'list' not callable error makes no sense anyway
-- cut and paste the same code to the prompt and it runs.

I notice it's important not to *close* the plt.plot window, if you want
to use it again same session.  Just change domain and range and replot,
without closing the plot window in between.

Here's my result:
http://www.inetarena.com/~pdx4d/ocn/python/cospoly.png  Notice this
polynomial crosses the x-axis 10 times!  But it's a degree 40 equation no?
No.  I've got 20 non-zero terms, but that starts from x^0, so I only get
up to x^38.  38th degree equation, with 38 roots.  A quick check in J
shows it has these 10 real roots (where it crosses x), plus 14 conjugate
pairs of complex roots (J has a built-in for finding these roots).

The real roots are:

         _14.1262           14.1262
         _10.9956           10.9956
          7.85398          _7.85398
          4.71239          _4.71239
           1.5708           _1.5708

which corresponds to plus/minus pi/2 3pi/2 5pi/2 7pi/2 and 9pi/2 --
except these values don't match the pi values.  Because this 38 degree
poly is an *approximation* (gotta have infinite terms to be spot on).

The only Python functions I added, to create a domain and range for
plt.plot(domain,range), were:

  def mkdomain(x):  # kinda ugly, but it works
      x *= 10
      return map(lambda x: x*0.1, range(-x,0,1)) + map(lambda x: x*0.1, 
range(x+1))

  def mkrange(d):
      return [cospoly(i,20) for i in d]  # the 20 fixes the number of terms

So this graph was a result of going:

   >>> from scipy import plt
   >>> import plotplay # the generator, cospoly, and the above two 
functions [1]
   >>> d = plotplay.mkdomain(16)  # from -16 to 16 in 0.1 increments
   >>> r = plotplay.mkrange (d)   # feed this domain to cospoly
   >>> plt.plot(d,r)

and that's it!

Kirby

[1] recall that we need to import some stuff in plotplay:

from __future__ import generators, division
from operator import add