ส่งงาน week11

week11/integrat.py

+#Trapezoidal rule
+# โจทย์ f(x) = e^x*sin(x) ,n=10
+from math import e,sin
+def f(x) :
+    return e**x*sin(x)
+#n คือ จำนวนกรแบ่งช่วง
+#h คือ ค่าของช่วงที่แบ่ง from (b-a)/n
+def trap(f,a,h,n) :
+    return sum([ (1/2)*h*(f(a+i*h)+f(a+(i+1)*h)) for i in range(n) ])
+def area(f,a,b,n):
+    return sum([ (1/2)*((b-a)/n)*(f(a+i*((b-a)/n))+f(a+(i+1)*((b-a)/n))) for i in range(n) ])
+trap(f,0.1,0.05,10),area(f,0.1,0.6,10)
+
+.............................................................................................
+#หา integrate จาก Sympy
+x,y=sympy.Symbol('x'),sympy.Symbol('y')
+y=3*x**2+9*x+2
+z=sympy.integrate(y,(x,0.1,0.6))
+from math import e
+zz = e**-x
+zzz = sympy.integrate(zz,(x,0,1))
+y,z,zzz
+
+#Simpson's rule
+
+import sympy
+from math import e
+f =e**(-x)
+sympy.integrate(f,(x,0,1))
+
+#print 0.632120558828558
+
+def simpson(f, a, b, n):
+    h = (b-a)/n
+    x0 = a # x_i = x_0 + i*h
+    #f(x0) + 4*f(x1) + 2*f(x2) + .... + 4*f(x_n-1) + f(x_n)
+    s = f(x0)
+    
+    # odd - x1, x3, x5, ... 
+    #s += 4*f(x1) + 4*f(x3) + 4*f(x5) + ...
+    for i in range(1,n,2):
+        s += 4*f(x0+i*h)
+        
+    # even - x2, x4, x6, ...
+    #s += 2*f(x2) + 2*f(x4) + 2*f(x6) + ...
+    for i in range(2,n-1,2):
+        s += 2*f(x0+i*h)
+    
+    s += f(x0+n*h)
+    return s*h/3
+
+from math import e
+from sympy import *
+def f(x): return e**(-x)
+a,b = 0,1
+r=0.632120558828558
+simpson(f, a, b, 4)
+for v in [ (n,simpson(f,a,b,n),area(f,a,b,n)) for n in range(4, 22, 2) ] :
+    print ('{:5} {:20} {:20}'.format(v[0],abs(v[1]-r),v[2]))