ModSim: Finished assignment 2.5

modsim/ass2/q5.c

@@ -42,7 +42,8 @@ double f2(double x) {
 }
 
 #define PRINT_INTEGRAL(func, method, a, b) (printf(#func " from " #a " to " \
-    #b " using %-19s %.11f\n", #method " method:", integral(&func, &method, a, b)))
+    #b " using %-19s %.11f\n", #method " method:", \
+            integral(&func, &method, a, b)))
 #define PRINT_GAUSS(func, a, b) (printf(#func " from " #a " to " \
     #b " using %-19s %.11f\n", "gauss method:", gauss(&func, a, b)))
 

modsim/ass2/q5.py

+# Gauss two-point integration formula from:
+# <http://kr.cs.ait.ac.th/~radok/math/mat7/step32.htm>
+
+from math import sqrt, exp, sin, pi
+
+def rectangle(fn, a, b):
+    return (b - a) * fn((a + b) / 2.0)
+
+def trapezoidal(fn, a, b):
+    return 0.5 * (b - a) * (fn(a) + fn(b))
+
+def simpson(fn, a, b):
+    return (2 * rectangle(fn, a, b) + trapezoidal(fn, a, b)) / 3.0
+
+def gauss(fn, a, b):
+    s = (b - a) / sqrt(3)
+    # calculate abscissae
+    left  = (b + a - s) / 2.0
+    right = (b + a + s) / 2.0
+    return (b - a) * (fn(left) + fn(right)) / 2
+
+def integrate(method, DX, fn, a, b):
+    surface = 0.0;
+    x = a
+
+    for i in xrange(int((b - a) / DX)):
+        x += DX
+        surface += method(fn, x, x + DX);
+
+    return surface;
+
+
+integrals = [(lambda x: exp(-x), 0, 1),
+        (lambda x: x * exp(-x), 0, 2),
+        (lambda x: x * exp(-x), 0, 20),
+        (lambda x: x * exp(-x), 0, 200),
+        (lambda x: sin(x), 0, 8 * pi)]
+
+methods = [rectangle, trapezoidal, simpson]
+
+DX = 1e-3
+
+for i in integrals:
+    for method in methods:
+        print 'integrate from %d to %d using %-11s: %.16f' \
+                % (i[1], i[2], method.func_name, integrate(method, DX, *i))
+    print 'integrate from %d to %d using %-11s: %.16f' \
+            % (i[1], i[2], 'gauss', gauss(*i))
+    print '-------------'
+