add: task 1, subtask 1-3 python reimplementation

Exercise_03/NSSC_1.m

@@ -1,28 +1,42 @@
-%% eulero forward DIrichlet
+%% Euler forward Dirichlet
-L=pi;
+L = 1;  % domain size (in the lecture notes this is denote h)
-T=2;
+T = 2;  % time limit (max time)
-f=@(x,t) 0*x.*t;
+f = @(x,t) 0*x.*t;    % rhs of the more general equation `u_t - d u_xx = f`
-c1=@(t) 1+0*t;
+c1 = @(t) 1+0*t;      % _right_ boundary condition
-c2=@(t) 0*t;
+c2 = @(t) 0*t;        % _left_ boundary condition
-u0=@(x) 0*x;
+u0 = @(x) 0*x;        % initial values
-D=0.5;
+D = 0.5;              % diffusion parameter `d` in `u_t - d u_xx = f`
-%uex=@(x,t) cos(x).*exp(t);
+%uex = @(x,t) cos(x).*exp(t);
 
-N=10;
+N = 10;               % nr. of _space_ discretization points
-K=200;
+K = 200;              % nr. of _time_ discretization points
-[x,t,u]=Dirichlet_EA(L,N,T,K,c1,c2,f,u0,D);
+[x, t, u] = Dirichlet_EA(L, N, T, K, c1, c2, f, u0, D);
 
-figure(1)
+% Report stability condition `D Delta T / (Delta x)^2 > 0.5`
-for ii=1:K+1 
+Delta_T = T / K;
-    plot(x,u(:,ii)');
+Delta_x = L / N;
-    xlim([0 L])
+d = D * Delta_T / Delta_x^2;
-    pause(0.05);
+fprintf("Stability Condition: 0.5 >= D * Delta_T / Delta_x^2 = %f\n", d)
+if d > 0.5
+    fprintf("-> NOT Stable\n")
+else
+    fprintf("-> Stable\n")
 end
 
-space=linspace(0,L,101);
+figure(1)
-time=linspace(0,T,201);
+for ii = 1:K+1 % iterates time
-[xx,yy]=meshgrid(time,space);
+    hold on
-%exsol=uex(yy,xx);
+    plot(x, u(:, ii)');
+    xlim([0 L])
+    pause(0.05);
+    hold off
+end
+
+% 3D plot of space solution over time
+space = linspace(0,L,101);
+time = linspace(0,T,201);
+[xx,yy] = meshgrid(time,space);
+%exsol = uex(yy,xx);
 figure(2)
 mesh(t,x,u)
 %figure(2)

Exercise_03/task01_1-2.py

@@ -0,0 +1,63 @@
+# Task 1.1, 1.2
+import numpy as np
+from typing import Callable
+from matplotlib import pyplot as plt
+
+# Config
+D = 1e-6            # diffusion coefficient
+h = 1               # space domain (max x size)
+T = 2e6             # solution end time
+nx = 50             # nr of space discretization points
+nt = 20000          # nr of time discretization points
+
+# derived constants
+dx = h / (nx - 1)   # space step size
+dt = T / (nt - 1)   # time step size
+d = dt * D / dx**2  # stability/stepsize coefficient
+
+# report stability
+if d > 0.5:
+    print("NOT Stable")
+else:
+    print("Stable")
+
+# explicit scheme integration for `u_t - D u_xx = 0` with boundary conditions
+# enforced by `set_bounds` and initial conditions `initial`.
+def integrate(*, name: str, initial: np.array, set_bounds: Callable[[np.array], None]) -> None:
+    C = initial
+    # Setup boundary conditions
+    set_bounds(C)
+
+    i = 0       # index for plot generation
+    plt.figure(figsize = (8, 6), dpi = 100)
+    for t in range(nt):
+        # every 400'th time step save a plot
+        if t % (nt // 400) == 0:
+            plt.clf()
+            plt.plot(np.linspace(0, h, nx), C)
+            plt.xlim([0, h])
+            plt.ylim([0, 1.2])
+            plt.savefig(f"plots/{name}_{i:0>5}.png")
+            i += 1
+        # update solution using the explicit schema
+        C[1:-1] += d * (C[2:] - 2 * C[1:-1] + C[:-2])
+        # update right Neumann BC
+        set_bounds(C)
+
+# Subtask 1 boundary conditions (Dirichlet and Neumann)
+def bounds_1(C):
+    C[0] = 1
+    C[-1] = C[-2]
+
+# Subtask 2 boundary conditions (two Dirichlet)
+def bounds_2(C):
+    C[0] = 1
+    C[-1] = 0
+
+# run simulations
+integrate(name = 'task01_1', initial = np.zeros(nx), set_bounds = bounds_1)
+integrate(name = 'task01_2', initial = np.zeros(nx), set_bounds = bounds_2)
+
+# to convert generated image sequence to video use:
+# $> ffmpeg -r 60 -i plots/task01_1_%05d.png -pix_fmt yuv420p video_1_1.mp4
+# $> ffmpeg -r 60 -i plots/task01_2_%05d.png -pix_fmt yuv420p video_1_2.mp4

Exercise_03/task01_3.py

@@ -0,0 +1,50 @@
+# Task 1.3
+import numpy as np
+from typing import Callable
+from matplotlib import pyplot as plt
+
+# Config
+D = 1e-6            # diffusion coefficient
+h = 1               # space domain (max x size)
+T = 2e6             # solution end time
+nx = 50             # nr of space discretization points
+nt = 20000          # nr of time discretization points
+
+# derived constants
+dx = h / (nx - 1)   # space step size
+dt = T / (nt - 1)   # time step size
+d = dt * D / dx**2  # stability/stepsize coefficient
+
+# Setup implicit scheme equation matrix
+T = (1 + 2 * d) * np.eye(nx) - d * np.eye(nx, k = 1) - d * np.eye(nx, k = -1)
+# fix boundary condition equations
+T[0, 0] = 1         # Left Dirichlet BC
+T[0, 1] = 0
+T[-1, -2] = 1       # Right Neumann BC
+T[-1, -1] = 0
+
+# Set initial solution
+C = np.zeros(nx)
+C[0] = 1
+C[-1] = C[-2] # (0 = 0)
+
+i = 0       # index for plot generation
+plt.figure(figsize = (8, 6), dpi = 100)
+for t in range(nt):
+    # every 400'th time step save a plot
+    if t % (nt // 400) == 0:
+        plt.clf()
+        plt.plot(np.linspace(0, h, nx), C)
+        plt.xlim([0, h])
+        plt.ylim([0, 1.2])
+        plt.savefig(f"plots/task01_3_{i:0>5}.png")
+        i += 1
+    # update solution using the explicit schema
+    C = np.linalg.solve(T, C)
+    # fix BC conditions (theoretically, they are set by the update but for
+    # stability reasons (numerical) we enforce the correct values)
+    C[0] = 1
+    C[-1] = C[-2]
+
+# to convert generated image sequence to video use:
+# $> ffmpeg -r 60 -i plots/task01_3_%05d.png -pix_fmt yuv420p video_1_3.mp4