99
1010def f (x ):
1111 """ the function we are integrating """
12- return 1.0 + x * 0.25 * np .sin (np .pi * x )
12+ #return 1.0 + x*0.25*np.sin(np.pi*x)
13+ return (1.0 + x * 0.25 * np .sin (5 * np .pi * x )) * np .exp (- (x - 1.0 )** 2 / 0.5 )
1314
1415
1516plt .rcParams .update ({'xtick.labelsize' : 18 ,
@@ -56,22 +57,31 @@ def rectangle(xp, fp, a, b):
5657
5758 ax = plt .gca ()
5859
60+ integral = 0.0
61+
5962 for n in range (len (xp )- 1 ):
6063
6164 xl = xp [n ]
6265 xr = xp [n + 1 ]
6366
64- print (xl )
6567 fl = fp [n ]
6668
6769 # shade region
6870 verts = [(xl , 0 ), (xl , fl ), (xr , fl ), (xr , 0 )]
6971 ax .add_patch (Polygon (verts , facecolor = "0.8" , edgecolor = "k" ))
7072
73+ # and bonus! actually compute the integral in this approximation
74+ integral += (xr - xl ) * fl
75+
76+ return integral
77+
78+
7179def trapezoid (xp , fp , a , b ):
7280
7381 ax = plt .gca ()
7482
83+ integral = 0.0
84+
7585 for n in range (len (xp )- 1 ):
7686
7787 xl = xp [n ]
@@ -84,10 +94,16 @@ def trapezoid(xp, fp, a, b):
8494 verts = [(xl , 0 ), (xl , fl ), (xr , fr ), (xr , 0 )]
8595 ax .add_patch (Polygon (verts , facecolor = "0.8" , edgecolor = "k" ))
8696
97+ integral += 0.5 * (xr - xl ) * (fl + fr )
98+
99+ return integral
100+
87101def simpsons (xp , fp , a , b ):
88102
89103 ax = plt .gca ()
90104
105+ integral = 0.0
106+
91107 for n in range (0 , len (xp )- 1 , 2 ):
92108
93109 # we need to handle the 1 bin case specially
@@ -124,13 +140,16 @@ def simpsons(xp, fp, a, b):
124140 verts = [(xl , 0 )] + simpvert + [(xr , 0 )]
125141 ax .add_patch (Polygon (verts , facecolor = "0.8" , edgecolor = "k" ))
126142
143+ integral += (xr - xl ) / 6.0 * (f0 + 4 * f1 + f2 )
144+
145+ return integral
127146
128147def main ():
129148
130149 N_fine = 200
131150
132151 # the number of bins to divide [a, b]
133- N_bins = 1
152+ N_bins = 16
134153
135154 xmin = 0.5
136155 xmax = 1.5
@@ -158,18 +177,18 @@ def main():
158177 plt .clf ()
159178
160179 plot_base (xp , fp , xfine , a , b )
161- rectangle (xp , fp , a , b )
180+ I_r = rectangle (xp , fp , a , b )
162181
163- plt .savefig ("rectangle.pdf " , bbox_inches = "tight" )
182+ plt .savefig (f"rectangle_N { N_bins } .png " , bbox_inches = "tight" )
164183
165184 # trapezoid method
166185
167186 plt .clf ()
168187
169188 plot_base (xp , fp , xfine , a , b )
170- trapezoid (xp , fp , a , b )
189+ I_t = trapezoid (xp , fp , a , b )
171190
172- plt .savefig ("trapezoid.pdf " , bbox_inches = "tight" )
191+ plt .savefig (f"trapezoid_N { N_bins } .png " , bbox_inches = "tight" )
173192
174193
175194 # simpsons method
@@ -190,11 +209,11 @@ def main():
190209 label_mid = label_mid )
191210
192211
193- simpsons (xp , fp , a , b )
194-
195- plt .savefig ("simpsons.pdf" , bbox_inches = "tight" )
212+ I_s = simpsons (xp , fp , a , b )
196213
214+ plt .savefig (f"simpsons_N{ N_bins } .png" , bbox_inches = "tight" )
197215
216+ print (f"integral approximations: { I_r } , { I_t } , { I_s } " )
198217
199218if __name__ == "__main__" :
200219 main ()
0 commit comments