Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import numpy as np
- import matplotlib as mpl
- import JUtil
- plot = JUtil.qjplot
- mpl.use("TKAgg")
- import matplotlib.pylab as plt
- plt.ion()
- def midpoint(f, x0, xn, n):
- X=[0.0]*(4*n)
- Y=[0.0]*(4*n)
- dx = float(xn-x0) / float(n)
- VLines = [x0+i*dx for i in xrange(n)]
- for i in xrange(n):
- X[4*i] = X[4*i+1] = float(x0+i*dx)
- X[4*i+2] = X[4*i+3] = x0+(i+1)*dx
- Y[4*i] = Y[4*i+3] = 0.0
- Y[4*i+1] = Y[4*i+2] = f( (X[4*i]+X[4*i+3]) / 2.0)
- return X,Y,VLines
- def trapezoid(f, x0, xn, n):
- X=[0.0]*(4*n)
- Y=[0.0]*(4*n)
- dx = float(xn-x0) / float(n)
- VLines = [x0+i*dx for i in xrange(n)]
- for i in xrange(n):
- X[4*i] = X[4*i+1] = float(x0+i*dx)
- X[4*i+2] = X[4*i+3] = x0+(i+1)*dx
- Y[4*i] = Y[4*i+3] = 0.0
- Y[4*i+1] = f(X[4*i])
- Y[4*i+2] = f(X[4*i+3])
- return X,Y,VLines
- def simpson(f, x0, xn, n):
- X,Ym,V = midpoint(f,x0,xn,n)
- Xt,Yt,Vt = trapezoid(f,x0,xn,n)
- Y = [(2.0*M+T)/3.0 for M,T in zip(Ym,Yt)]
- return X,Y,V
- if '__main__' in __name__:
- a=0.0
- b=1.0
- x = lambda n: np.linspace(a,b,n)
- f = lambda x: np.cos(np.pi*x)+x
- plt.clf()
- #Origin
- plt.hlines(0.0,-10,10, 'k')
- plt.vlines(0.0, -10, 10, 'k')
- #Midpoint
- X,Y,VLines = midpoint(f,a,b,4)
- plot(X,Y, 'r', rebound=0)
- #Trapezoid
- X,Y,VLines = trapezoid(f,a,b,4)
- plot(X,Y, 'g', rebound=0)
- #Simpson
- X,Y,VLines = simpson(f,a,b,4)
- plot(X,Y, 'b', rebound=0)
- #Grid lines
- plt.vlines(VLines, 0, max(f(x(200))), 'k', linestyles='dashed')
- #Plot function
- plot(x(200),f(x(200)), 'k', rebound=1)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement