Hello all,

I've been working on writing some PyX code in order to
illustrate the different methods of numerical
integration. So far I've succeeded in writing code for
all of the different Riemann sums, and the trapezoidal
rule. For me this just leaves Simpson's rule to work
on (which unfortunately is quite a bit more
complicated). Here is my work in progress code:

from math import *
from pyx import *
from pyx.graph import axis

#Define Variables
def f(x):
    return sin(x) + x

x0,y0 = 0,0     # Starting point for Riemann sum
intl = 4        # width of total interval
nsubintl = 6    # number of subintervals

# Code

w0 = intl/nsubintl # width of subinterval

g = graph.graphxy(width=8,
        x=axis.linear(min=0, max=intl, title="$x$"),
        y=axis.linear(title="$\sin(x) + x$"))

g.plot(graph.data.function("y(x)=f(x)",
context=locals(), points=500))

for i in range(0, nsubintl):
    def g(x):
        if x0+i*w0 <= x <= x0+(i+1)*w0: return
(f(x0+i*w0) + f(x0+(i+0.5)*w0) + f(x0+(i+1)*w0) )/(2 *
w0 * w0) * x * x + (f(x0+(i+1)*w0) - f(x0+i*w0)  )/(2
* w0) * x + f(x0+(i+0.5)*w0)

    g.plot(graph.data.function("y(x)=g(x)",
context=locals(), points=500))

else:
    g.dolayout()
    g.stroke(g.ygridpath(0),
[style.linestyle.dashed,style.linewidth.Thin])
    g.stroke(g.xgridpath(0),
[style.linestyle.dashed,style.linewidth.Thin])
    g.writeEPSfile("trapezoidal_rule")

The problem is that when I run it I get the error:

bsd% python simpsons_rule.py 
Traceback (most recent call last):
  File "simpsons_rule.py", line 27, in <module>
    g.plot(graph.data.function("y(x)=g(x)",
context=locals(), points=500))
AttributeError: 'function' object has no attribute
'plot'

Can anybody shed some light on the problem?

As long as I'm on the subject, could anybody help me
actually figure out how to shade the regions required
for Simpson's rule? In each region I will want to
connect three straight lines, with a parabola on top,
and then fill the total area with a solid color. When
I coded the trapezoidal rule, I used:

    trap = path.path(path.moveto(x1+i*w1,y1),
path.lineto(x1+i*w1, b),
                     path.lineto(x1+(i+1)*w1,c),
                     path.lineto(x1+(i+1)*w1,y1),
                     path.closepath())
    g.stroke(trap,[deco.filled([color.gray(0.8)]) ])

but I'm not sure how to actually include a curved
portion (the parabola) as part of a path. Any
suggestions would be greatly appreciated.

For reference here is the working code I'm using for
the left Riemann sum:

from math import *
from pyx import *
from pyx.graph import axis

#Define Variables
def f(x):
    return sin(x) + x

x0,y0 = 0,0     # Starting point for Riemann sum
intl = 2*pi     # width of total interval
nsubintl = 20   # number of subintervals

# Code

g = graph.graphxy(width=8,
        x=axis.linear(min=0, max=intl, title="$x$"),
        y=axis.linear(title="$\sin(x) + x$"))

g.plot(graph.data.function("y(x)=f(x)",
context=locals(), points=500))

g.dolayout()

g.stroke(g.ygridpath(0),
[style.linestyle.dashed,style.linewidth.Thin])
g.stroke(g.xgridpath(0),
[style.linestyle.dashed,style.linewidth.Thin])

# Left Riemann Sum

x1,y1 = g.pos(x0,y0) # Starting point on the graph
itself
x2,y2 = g.pos(x0+intl/nsubintl, 0) # Need to find out
the distance of intl/subintl on the graph

w0 = intl/nsubintl # width of subinterval
w1= x2-x1 # width of subinterval on the graph 

# rect2 = path.rect(a, b, c, d), (a,b) origin, c
width, d height.

for i in range(0, nsubintl):
    a,b = g.pos(x0,f(x1+i*w0)) # Get the height of the
function b (a is useless)
    h = b - y1
   
g.stroke(path.rect(x1+i*w1,y1,w1,h),[deco.filled([color.gray(0.8)])
])
else:
    g.writeEPSfile("left_riemann_sum")


It's probably not the neatest or best way to achieve
my goals, but it works nonetheless. 

Many thanks,
-Daniel


      

-------------------------------------------------------------------------
This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
Build the coolest Linux based applications with Moblin SDK & win great prizes
Grand prize is a trip for two to an Open Source event anywhere in the world
http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________
PyX-user mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/pyx-user

Reply via email to