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