Advertisement
Guest User

Numerical Analysis

a guest
Dec 8th, 2010
122
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 1.59 KB | None | 0 0
  1. import numpy as np
  2. import matplotlib as mpl
  3. import JUtil
  4. plot = JUtil.qjplot
  5. mpl.use("TKAgg")
  6. import matplotlib.pylab as plt
  7. plt.ion()
  8.  
  9. def midpoint(f, x0, xn, n):
  10.     X=[0.0]*(4*n)
  11.     Y=[0.0]*(4*n)
  12.    
  13.     dx = float(xn-x0) / float(n)
  14.     VLines = [x0+i*dx for i in xrange(n)]
  15.     for i in xrange(n):
  16.         X[4*i] = X[4*i+1] = float(x0+i*dx)
  17.         X[4*i+2] = X[4*i+3] = x0+(i+1)*dx
  18.  
  19.         Y[4*i] = Y[4*i+3] = 0.0
  20.         Y[4*i+1] = Y[4*i+2] = f( (X[4*i]+X[4*i+3]) / 2.0)
  21.     return X,Y,VLines
  22.  
  23. def trapezoid(f, x0, xn, n):
  24.     X=[0.0]*(4*n)
  25.     Y=[0.0]*(4*n)
  26.    
  27.     dx = float(xn-x0) / float(n)
  28.     VLines = [x0+i*dx for i in xrange(n)]
  29.     for i in xrange(n):
  30.         X[4*i] = X[4*i+1] = float(x0+i*dx)
  31.         X[4*i+2] = X[4*i+3] = x0+(i+1)*dx
  32.  
  33.         Y[4*i] = Y[4*i+3] = 0.0
  34.         Y[4*i+1] = f(X[4*i])
  35.         Y[4*i+2] = f(X[4*i+3])
  36.     return X,Y,VLines
  37.  
  38. def simpson(f, x0, xn, n):
  39.     X,Ym,V = midpoint(f,x0,xn,n)
  40.     Xt,Yt,Vt = trapezoid(f,x0,xn,n)
  41.     Y = [(2.0*M+T)/3.0 for M,T in zip(Ym,Yt)]
  42.     return X,Y,V
  43.  
  44. if '__main__' in __name__:
  45.     a=0.0
  46.     b=1.0
  47.     x = lambda n: np.linspace(a,b,n)
  48.     f = lambda x: np.cos(np.pi*x)+x
  49.  
  50.     plt.clf()
  51.     #Origin
  52.     plt.hlines(0.0,-10,10, 'k')
  53.     plt.vlines(0.0, -10, 10, 'k')
  54.  
  55.     #Midpoint
  56.     X,Y,VLines = midpoint(f,a,b,4)
  57.     plot(X,Y, 'r', rebound=0)
  58.  
  59.     #Trapezoid
  60.     X,Y,VLines = trapezoid(f,a,b,4)
  61.     plot(X,Y, 'g', rebound=0)
  62.  
  63.     #Simpson
  64.     X,Y,VLines = simpson(f,a,b,4)
  65.     plot(X,Y, 'b', rebound=0)
  66.  
  67.     #Grid lines
  68.     plt.vlines(VLines, 0, max(f(x(200))), 'k', linestyles='dashed')
  69.  
  70.     #Plot function
  71.     plot(x(200),f(x(200)), 'k', rebound=1)
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement