Finished Exercise 1

week1/ex1.py

@@ -1,36 +1,125 @@
 #!/usr/bin/env python3
 
-import numpy as np
-import matplotlib.pyplot as plt
-
 
 # Truncation of Euler's Constant
-e_OEIS = 2.7182818284590452353
-N_max = 25
-N_range = np.arange(0, N_max)
 
-# a) Simple Relative Error Determine
+def mainA():
-e_approx = np.zeros(N_max)
+    e_OEIS = 2.7182818284590452353
-rel_err = np.zeros(N_max)
+    N_max = 25
+    N_range = np.arange(0, N_max)
     
-prev_e = 0
+    # a) Simple Relative Error Determination
-for i, n in enumerate(N_range):
+    e_approx = np.empty(N_max)
-    e_approx[i] = prev_e + 1/np.math.factorial(n)
+    rel_err = np.empty(N_max)
     
-    rel_err[i] = abs( e_approx[i] - e_OEIS ) / e_OEIS
+    prev_e = 0
-    prev_e = e_approx[i]
+    for i, n in enumerate(N_range):
+        e_approx[i] = prev_e + 1/np.math.factorial(n)
     
-print(e_approx)
+        rel_err[i] = abs( e_approx[i] - e_OEIS ) / e_OEIS
-
+        prev_e = e_approx[i]
-fig, ax = plt.subplots()
-ax.plot(N_range, rel_err)
-ax.set_yscale('log')
-ax.grid()
-ax.set_ylabel("Relative Error")
-ax.set_xlabel("N")
-
-plt.show()
     
 
-# b) Varying Floating Point Precision
+    ## Plot it all
+    fig, ax = plt.subplots()
+    ax.set_title("Relative Error for std float")
+    ax.plot(N_range, rel_err)
+    ax.set_yscale('log')
+    ax.grid()
+    ax.set_ylabel("Relative Error")
+    ax.set_xlabel("N")
     
+    plt.show()
+    
+def mainB():
+    # b) Varying Floating Point Precision
+    e_OEIS = 2.7182818284590452353
+    N_max = 25
+    N_range = np.arange(0, N_max)
+
+    e_approx = np.zeros(N_max)
+    rel_err = np.zeros(N_max)
+    
+    series_element_64 = np.empty(N_max, dtype=np.float64)
+    series_element_32 = np.empty(N_max, dtype=np.float32)
+    
+    series64_diff = np.empty(N_max, dtype=np.float64)
+    series32_diff = np.empty(N_max, dtype=np.float64)
+    
+    ## Determine the coefficients
+    for i,n in enumerate(N_range):
+        e_n = 1/np.math.factorial(n)
+    
+        series_element_64[i] = e_n # auto casts to float64
+        series_element_32[i] = e_n # auto casts to float32
+    
+    ## Make a Cumulative Sum of the elements
+    series64 = np.cumsum(series_element_64)
+    series32 = np.cumsum(series_element_32)
+    
+    ## Determine the relative errror
+    for i,n in enumerate(N_range):
+        series64_diff[i] = abs( series64[i] - e_OEIS ) / e_OEIS
+        series32_diff[i] = abs( series32[i] - e_OEIS ) / e_OEIS
+    
+
+    ## Plot it all
+    fig, ax = plt.subplots()
+    ax.set_title("Relative Error for float64 and float32")
+    ax.plot(N_range, series64_diff, label="float64")
+    ax.plot(N_range, series32_diff, label="float32")
+    ax.set_yscale('log')
+    ax.grid()
+    ax.legend()
+    ax.set_ylabel("Relative Error")
+    ax.set_xlabel("N")
+    
+    plt.show()
+    
+    
+def mainC():
+    # c) Relative errors for different rounding accuracies
+    round_accuracies = np.arange(1,6)
+
+    e_OEIS = 2.7182818284590452353
+    N_max = 25
+    N_range = np.arange(0, N_max)
+    
+    # a) Simple Relative Error Determination
+    e_approx = np.empty(N_max)
+    
+    rel_err = np.empty( (len(round_accuracies), N_max) )
+
+    for i,d in enumerate(round_accuracies):
+        prev_e = 0
+        for j, n in enumerate(N_range):
+            e_approx[i] = prev_e + round(1/np.math.factorial(n), d)
+    
+            rel_err[i,j] = abs( e_approx[i] - e_OEIS ) / e_OEIS
+            prev_e = e_approx[i]
+    
+
+    ## Plot it all
+    fig, ax = plt.subplots()
+    ax.set_title("Relative Errors for std float with rounding at various digits d")
+    for i,d in enumerate(round_accuracies):
+        ax.plot(N_range, rel_err[i], label="d = {}".format(d))
+
+    ax.set_yscale('log')
+    ax.grid()
+    ax.set_ylabel("Relative Error")
+    ax.set_xlabel("N")
+    ax.legend()
+    
+    plt.show()
+
+
+
+
+if __name__ == "__main__":
+    import numpy as np
+    import matplotlib.pyplot as plt
+
+    mainA()
+    mainB()
+    mainC()