add ex2b: numerical integration

Exercise2b/numerical_integration.py

@@ -0,0 +1,56 @@
+# Exercise 2b
+
+# Here we have a Python function which calculates the integral in two
+# different ways: numerically and analytically.
+#
+# We want to check the precision of the numerical integration as a function of
+# the number of steps (subintervals). To do this, we calculate and print the
+# relative differences between the analytic result and the numerical result
+# for different values of the number of steps.
+#
+# Steps:
+# 0. Familizare yourselves with code below.
+# 1. Implement the serial version using a for loop or map
+# 2. Implement the parallel version using multiprocessing.Pool
+# 3. Time both versions
+# 4. What (if any) do you get?
+
+def integrate(f, a, b, n):
+    "Perform numerical integration of f in range [a, b], with n steps"
+    s = []
+    for i in range(n):
+        dx = (b - a) / n
+        x = a + (i + 0.5) * dx
+        y = f(x)
+        s = s + [y * dx]
+    return sum(s)
+
+def f(x):
+    "A polynomial that we'll integrate"
+    return x ** 4 - 3 * x
+
+def F(x):
+    "The analatic integral of f. (F' = f)"
+    return 1 / 5 * x ** 5 - 3 / 2 * x ** 2
+
+def compute_error(n):
+    "Calculate the difference between the numerical and analytical integration results"
+    a = -1.0
+    b = +2.0
+    F_analytical = F(b) - F(a)
+    F_numerical = integrate(f, a, b, n)
+    return abs((F_numerical - F_analytical) / F_analytical)
+
+def main():
+    ns = [10_000, 25_000, 50_000, 75_000]
+    errors = ... # TODO: write a for loop, serial map, and parallel map here
+
+    for n, e in zip(ns, errors):
+        print(f'{n} {e:.8%}')
+
+if __name__ == '__main__':
+    main()
+
+# Bonus steps, very optional:
+# 6. Implement a parallel version with threads (using multiprocessing.pool.ThreadPool).
+# 7. Time this version, and hypothetize about the result.

Exercise2b/numerical_integration_solution.py

@@ -0,0 +1,30 @@
+import sys
+from numerical_integration import compute_error
+
+def main(arg):
+    ns = [10_000, 25_000, 50_000, 75_000]
+
+    match arg:
+        case 'for':
+            errors = []
+            for n in ns:
+                errors += [compute_error(n)]
+        case 'lc':
+            errors = [compute_error(n) for n in ns]
+        case 'map':
+            errors = list(map(compute_error, ns))
+        case 'mp':
+            from multiprocessing import Pool as ProcessPool
+            with ProcessPool() as pool:
+                errors = pool.map(compute_error, ns)
+        case 'mt':
+            from multiprocessing.pool import ThreadPool
+            with ThreadPool(10) as pool:
+                errors = pool.map(compute_error, ns)
+
+    for n, e in zip(ns, errors):
+        print(f'{n} {e:.8%}')
+
+if __name__ == '__main__':
+    arg = (sys.argv[1:] + ['for'])[0]
+    main(arg)