From aa99b82a193c696460ac63ad77d5ca7dc91ea3ec Mon Sep 17 00:00:00 2001
From: David Doebel <ddoebel@ethz.ch>
Date: Tue, 3 Mar 2026 21:55:01 +0100
Subject: [PATCH] Add exercise NumFin Series 1

---
 .../NumFinSeries1/1_exercise3_template.py     | 109 ++++++++++++++++++
 main.py                                       |   0
 2 files changed, 109 insertions(+)
 create mode 100644 NumFinSeries/NumFinSeries1/1_exercise3_template.py
 create mode 100644 main.py

diff --git a/NumFinSeries/NumFinSeries1/1_exercise3_template.py b/NumFinSeries/NumFinSeries1/1_exercise3_template.py
new file mode 100644
index 0000000..9c55a73
--- /dev/null
+++ b/NumFinSeries/NumFinSeries1/1_exercise3_template.py
@@ -0,0 +1,109 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import numpy.linalg as lin
+
+
+# Set how floating-point errors are handled.
+np.seterr(all='raise')
+
+
+def initial_value(x):
+    return np.sin(np.pi/2 * x)
+
+
+#### exact solution at t=1 ####
+def exact_solution_at_1(x):
+
+    return np.exp(- np.pi ** 2 /4 * np.sin(np.pi / 2 * x))
+
+
+#### numerical scheme ####
+def eulerexplicit(N, M):
+
+    # Initialize the array u
+    u = np.linspace(0, 1, N)
+    u = initial_value(u)
+    # Initialize the G matrix
+    G = np.diag(np.repeat(-2,N),0) + np.diag(np.ones(N-1),-1) \
+        + np.diag(np.ones(N-1),1)
+    G[-1,-2] = 2
+    print(G)
+    # Initialize nu
+    k = 1/M
+    h = 1/N
+    nu = k/h**2
+    # Compute the solution
+    u = (np.identity(N) + nu * G) @ u
+    return u
+
+
+def eulerimplicit(N, M):
+
+    # Initialize the array u
+    u = np.linspace(0, 1, N)
+    u = initial_value(u)
+    # Initialize the G matrix
+    G = np.diag(np.repeat(-2, N), 0) + np.diag(np.ones(N - 1), -1) \
+        + np.diag(np.ones(N - 1), 1)
+    G[-1, -2] = 2
+    # Initialize nu
+    k = 1 / M
+    h = 1 / N
+    nu = k / h ** 2
+    # Compute the solution
+    u = np.linalg.solve(np.identity(N) - nu * G, u)
+    return u
+
+
+#### error analysis ####
+nb_samples = 5
+N = np.pow(2,np.arange(2,7))
+M = 2 * np.pow(4,np.arange(2,7))
+print(N)
+l2errorexplicit = np.zeros(nb_samples)  # error vector for explicit method
+l2errorimplicit = np.zeros(nb_samples)  # error vector for implicit method
+h2k = 1 / (N ** 2) + 1 / M
+
+
+#### Do not change any code below! ####
+try:
+    for i in range(nb_samples):
+        l2errorexplicit[i] = (1 / N[i]) ** (1 / 2) * lin.norm(
+            exact_solution_at_1(np.linspace(0, 1, N[i] + 1)[1:]) - eulerexplicit(N[i], M[i]), ord=2)
+    conv_rate = np.polyfit(np.log(h2k), np.log(l2errorexplicit), deg=1)
+    if np.isnan(conv_rate[0]):
+        raise Exception("Error unbounded for explicit method. Plots not shown.")
+    print("Explicit method converges: Convergence rate in discrete $L^2$ norm with respect to $h^2+k$: " + str(
+        conv_rate[0]))
+    plt.figure(figsize=[10, 6])
+    plt.loglog(h2k, l2errorexplicit, '-x', label='error')
+    plt.loglog(h2k, h2k, '--', label='$O(h^2+k)$')
+    plt.title('$L^2$ convergence rate for explicit method', fontsize=13)
+    plt.xlabel('$h^2+k$', fontsize=13)
+    plt.ylabel('error', fontsize=13)
+    plt.legend()
+    plt.plot()
+except Exception as e:
+    print(f"Exception: {e}")
+
+try:
+    for i in range(nb_samples):
+        l2errorimplicit[i] = (1 / N[i]) ** (1 / 2) * lin.norm(
+            exact_solution_at_1(np.linspace(0, 1, N[i] + 1)[1:]) - eulerimplicit(N[i], M[i]), ord=2)
+    conv_rate = np.polyfit(np.log(h2k), np.log(l2errorimplicit), deg=1)
+    if np.isnan(conv_rate[0]):
+        raise Exception("Error unbounded for implicit method. Plots not shown.")
+    print("Implicit method converges: Convergence rate in discrete $L^2$ norm with respect to $h^2+k$: " + str(
+        conv_rate[0]))
+    plt.figure(figsize=[10, 6])
+    plt.loglog(h2k, l2errorimplicit, '-x', label='error')
+    plt.loglog(h2k, h2k, '--', label='$O(h^2+k)$')
+    plt.title('$L^2$ convergence rate for implicit method', fontsize=13)
+    plt.xlabel('$h^2+k$', fontsize=13)
+    plt.ylabel('error', fontsize=13)
+    plt.legend()
+    plt.plot()
+except Exception as e:
+    print(f"Exception: {e}")
+
+plt.show()
\ No newline at end of file
diff --git a/main.py b/main.py
new file mode 100644
index 0000000..e69de29