init

src/consts.py

@@ -0,0 +1,5 @@
+l = 20
+T = 100
+
+D = 0.6
+H = 4

src/runnables/compare.py

@@ -0,0 +1,86 @@
+import itertools
+
+import util
+from solvers import solvers, solve_explicit, solve_implicit, solve_implicit_improved, solve_crank_nicolson, \
+  solve_crank_nicolson_improved
+from solvers.theoretical import solve_theoretical
+from util import Grid
+
+import numpy as np
+
+
+def test_order(errors, next_params):
+  errors = dict(errors)
+  for params in list(errors.keys()):
+    if params not in errors:
+      continue
+
+    # chain start
+    param_chain = []
+    error_chain = []
+
+    while params in errors:
+      param_chain.append(params)
+      error_chain.append(errors[params])
+
+      del errors[params]
+
+      params = next_params(params)
+
+    if len(param_chain) < 2:
+      continue
+
+    print(' ', ' -> '.join(map(str, param_chain)))
+    e = np.array(error_chain, dtype=float)
+    ed = e[:-1] / e[1:]
+    print(' ', ' '.join([f'{i:6.2f}' for i in ed]))
+
+
+def test_solver(solver, Is, Ks):
+  param_pairs = itertools.product(Is, Ks)
+  param_pairs = sorted(param_pairs, key=np.prod)
+
+  errors_mse = {}
+  errors_mae = {}
+  errors_mae_last = {}
+
+  for I, K in param_pairs:
+    grid = Grid(I, K)
+
+    solution = solver(grid)
+    if solution is None:
+      continue
+
+    reference = solve_theoretical(grid)
+
+    errors_mse[(I, K)] = util.mse(reference, solution)
+    errors_mae[(I, K)] = util.mae(reference, solution)
+    errors_mae_last[(I, K)] = util.mae_last(reference, solution)
+
+  print()
+  print()
+  print(solver.__name__)
+
+  print('O(h_x, h_t) MAEL       SHOULD BE >=2')
+  test_order(errors_mae_last, next_params=lambda params: (params[0] * 2, params[1] * 2))
+
+  print(f'O(h_x, h_t^2) MAEL    SHOULD BE >=4')
+  test_order(errors_mae_last, next_params=lambda params: (params[0] * 4, params[1] * 2))
+
+  print(f'O(h_x^2, h_t) MAEL    SHOULD BE >=4')
+  test_order(errors_mae_last, next_params=lambda params: (params[0] * 2, params[1] * 4))
+
+
+if __name__ == '__main__':
+  # Is = [4, 8, 16, 32, 64, 128, 256, 512]
+  # Ks = [4, 8, 16, 32, 64, 128, 256, 512]
+  # Is = [4, 8, 16, 32, 64, 128, 256, 512]
+  # Ks = [32, 64, 128, 256, 512, 1024, 2048, 4096]
+  Is = [4, 8, 16, 32, 64]
+  Ks = [32, 64, 128, 256, 512, 1024]
+
+  test_solver(solve_explicit, Is, Ks)
+  test_solver(solve_implicit, Is, Ks)
+  test_solver(solve_implicit_improved, Is, Ks)
+  test_solver(solve_crank_nicolson, Is, Ks)
+  test_solver(solve_crank_nicolson_improved, Is, Ks)

src/runnables/run_all.py

@@ -0,0 +1,20 @@
+import util
+from solvers import solvers
+from solvers.theoretical import solve_theoretical
+from util import Grid
+
+if __name__ == '__main__':
+  grid = Grid(16, 128)
+
+  reference = solve_theoretical(grid)
+
+  print()
+  print(grid)
+
+  for solver in solvers:
+    solution = solver(grid)
+
+    mse = util.mse(reference, solution)
+    mae = util.mae(reference, solution)
+
+    print(f'{solver.__name__:30} {mse:.4e}   {mae:.4e}')

src/runnables/where_mae.py

@@ -0,0 +1,27 @@
+from solvers import solvers, solve_theoretical
+from util import Grid
+
+import numpy as np
+
+if __name__ == '__main__':
+  I = 16
+  K = 512
+
+  grid = Grid(I, K)
+
+  for solver in solvers:
+    solution = solver(grid)
+    if solution is None:
+      continue
+
+    print(solver.__name__)
+
+    reference = solve_theoretical(grid)
+
+    e = np.abs(reference - solution)
+
+    flat_indices = np.argsort(e, axis=None)
+
+    for index in flat_indices[-10:]:
+      i, j = np.unravel_index(index, e.shape)
+      print(f'{i:4} {j:4} {e[i, j]}')

src/solvers/__init__.py

@@ -0,0 +1,12 @@
+from solvers.crank_nicolson import solve_crank_nicolson, solve_crank_nicolson_improved
+from solvers.explicit import solve_explicit
+from solvers.implicit import solve_implicit, solve_implicit_improved
+from solvers.theoretical import solve_theoretical
+
+solvers = [
+  solve_explicit,
+  solve_implicit,
+  solve_implicit_improved,
+  solve_crank_nicolson,
+  solve_crank_nicolson_improved,
+]

src/solvers/crank_nicolson.py

@@ -0,0 +1,102 @@
+import numpy as np
+import scipy.sparse
+import scipy.sparse.linalg
+
+import consts
+from util import Grid
+
+
+#@cachier(stale_after=datetime.timedelta(hours=1))
+def solve_crank_nicolson(grid: Grid) -> np.ndarray:
+  ug = np.zeros((grid.K + 1, grid.I + 1), dtype=float)
+
+  ug[0, :] = np.sin(np.pi * grid.xs / consts.l) ** 2 * 4 / 10
+
+  for k in range(1, grid.K + 1):
+    diags = {
+      -1: np.zeros(grid.I + 1, dtype=float),
+      0: np.zeros(grid.I + 1, dtype=float),
+      +1: np.zeros(grid.I + 1, dtype=float),
+    }
+    b = np.zeros(grid.I + 1, dtype=float)
+
+    # left
+    diags[0][0] = consts.H * grid.h_x + 1
+    diags[+1][0] = -1
+
+    # middle
+    for i in range(1, grid.I):
+      diags[-1][i] = grid.gamma
+      diags[0][i] = - 2 * grid.gamma - 2
+      diags[+1][i] = grid.gamma
+      b[i] = -(
+          ug[k - 1][i - 1] * grid.gamma +
+          ug[k - 1][i] * (-2 * grid.gamma + 2) +
+          ug[k - 1][i + 1] * grid.gamma
+      )
+
+    # right
+    diags[0][grid.I] = consts.H * grid.h_x + 1
+    diags[-1][grid.I] = -1
+
+    A = scipy.sparse.diags([
+      diags[-1][1:],
+      diags[0],
+      diags[+1][:-1],
+    ], [-1, 0, 1], format='csc')
+
+    ug[k, :] = scipy.sparse.linalg.spsolve(A, b)
+
+  return ug
+
+
+#@cachier(stale_after=datetime.timedelta(hours=1))
+def solve_crank_nicolson_improved(grid: Grid) -> np.ndarray:
+  ug = np.zeros((grid.K + 1, grid.I + 1), dtype=float)
+
+  ug[0, :] = np.sin(np.pi * grid.xs / consts.l) ** 2 * 4 / 10
+
+  for k in range(1, grid.K + 1):
+    diags = {
+      -1: np.zeros(grid.I + 1, dtype=float),
+      0: np.zeros(grid.I + 1, dtype=float),
+      +1: np.zeros(grid.I + 1, dtype=float),
+    }
+    b = np.zeros(grid.I + 1, dtype=float)
+
+    # left
+    diags[0][0] = consts.H * grid.gamma * grid.h_x + grid.gamma + 1
+    diags[+1][0] = - grid.gamma
+    b[0] = (
+        (- consts.H * grid.gamma * grid.h_x - grid.gamma + 1) * ug[k - 1, 0] +
+        grid.gamma * ug[k - 1, 1]
+    )
+
+    # middle
+    for i in range(1, grid.I):
+      diags[-1][i] = grid.gamma
+      diags[0][i] = - 2 * grid.gamma - 2
+      diags[+1][i] = grid.gamma
+      b[i] = -(
+          ug[k - 1][i - 1] * grid.gamma +
+          ug[k - 1][i] * (-2 * grid.gamma + 2) +
+          ug[k - 1][i + 1] * grid.gamma
+      )
+
+    # right
+    diags[0][grid.I] = consts.H * grid.gamma * grid.h_x + grid.gamma + 1
+    diags[-1][grid.I] = - grid.gamma
+    b[grid.I] = (
+        (- consts.H * grid.gamma * grid.h_x - grid.gamma + 1) * ug[k - 1, grid.I] +
+        grid.gamma * ug[k - 1, grid.I - 1]
+    )
+
+    A = scipy.sparse.diags([
+      diags[-1][1:],
+      diags[0],
+      diags[+1][:-1],
+    ], [-1, 0, 1], format='csc')
+
+    ug[k, :] = scipy.sparse.linalg.spsolve(A, b)
+
+  return ug

src/solvers/explicit.py

@@ -0,0 +1,32 @@
+from typing import Optional
+
+import numpy as np
+
+import consts
+from util import Grid
+
+
+#@cachier(stale_after=datetime.timedelta(hours=1))
+def solve_explicit(grid: Grid) -> Optional[np.ndarray]:
+  if grid.gamma > 1 / 2:
+    return None
+
+  ug = np.zeros((grid.K + 1, grid.I + 1), dtype=float)
+
+  ug[0, :] = np.sin(np.pi * grid.xs / consts.l) ** 2 * 4 / 10
+
+  for k in range(1, grid.K + 1):
+    # middle
+    ug[k, 1: grid.I] = (
+        grid.gamma * ug[k - 1, 0: grid.I - 1] +
+        (-2 * grid.gamma + 1) * ug[k - 1, 1: grid.I] +
+        grid.gamma * ug[k - 1, 2: grid.I + 1]
+    )
+
+    # left
+    ug[k, 0] = ug[k, 1] / (consts.H * grid.h_x + 1)
+
+    # right
+    ug[k, grid.I] = ug[k, grid.I - 1] / (consts.H * grid.h_x + 1)
+
+  return ug

src/solvers/implicit.py

@@ -0,0 +1,95 @@
+import numpy as np
+import scipy.sparse
+import scipy.sparse.linalg
+
+import consts
+
+
+#@cachier(stale_after=datetime.timedelta(hours=1))
+def solve_implicit(grid) -> np.ndarray:
+  ug = np.zeros((grid.K + 1, grid.I + 1), dtype=float)
+
+  ug[0, :] = np.sin(np.pi * grid.xs / consts.l) ** 2 * 4 / 10
+
+  for k in range(1, grid.K + 1):
+    diags = {
+      -1: np.zeros(grid.I + 1, dtype=float),
+      0: np.zeros(grid.I + 1, dtype=float),
+      +1: np.zeros(grid.I + 1, dtype=float),
+    }
+    b = np.zeros(grid.I + 1, dtype=float)
+
+    # left
+    diags[0][0] = consts.H * grid.h_x + 1
+    diags[+1][0] = -1
+
+    # middle
+    for i in range(1, grid.I):
+      diags[-1][i] = grid.gamma
+      diags[0][i] = - 2 * grid.gamma - 1
+      diags[+1][i] = grid.gamma
+      b[i] = -ug[k - 1, i]
+
+    # right
+    diags[-1][grid.I] = -1
+    diags[0][grid.I] = consts.H * grid.h_x + 1
+
+    A = scipy.sparse.diags([
+      diags[-1][1:],
+      diags[0],
+      diags[+1][:-1],
+    ], [-1, 0, 1], format='csc')
+
+    ug[k, :] = scipy.sparse.linalg.spsolve(A, b)
+
+  return ug
+
+
+#@cachier(stale_after=datetime.timedelta(hours=1))
+def solve_implicit_improved(grid) -> np.ndarray:
+  ug = np.zeros((grid.K + 1, grid.I + 1), dtype=float)
+
+  ug[0, :] = np.sin(np.pi * grid.xs / consts.l) ** 2 * 4 / 10
+
+  for k in range(1, grid.K + 1):
+    diags = {
+      -1: np.zeros(grid.I + 1, dtype=float),
+      0: np.zeros(grid.I + 1, dtype=float),
+      +1: np.zeros(grid.I + 1, dtype=float),
+    }
+    b = np.zeros(grid.I + 1, dtype=float)
+
+    # left
+    diags[0][0] = (
+        2 * consts.H * grid.gamma * grid.h_x +
+        2 * grid.gamma +
+        1
+    )
+    diags[+1][0] = - 2 * grid.gamma
+    b[0] = ug[k - 1, 0]
+
+    # middle
+    for i in range(1, grid.I):
+      diags[-1][i] = grid.gamma
+      diags[0][i] = - 2 * grid.gamma - 1
+      diags[+1][i] = grid.gamma
+      b[i] = -ug[k - 1, i]
+
+    # right
+    diags[0][grid.I] = (
+        2 * consts.H * grid.gamma * grid.h_x +
+        2 * grid.gamma +
+        1
+    )
+    diags[-1][grid.I] = - 2 * grid.gamma
+    b[grid.I] = ug[k - 1, grid.I]
+
+    A = scipy.sparse.diags([
+      diags[-1][1:],
+      diags[0],
+      diags[+1][:-1],
+    ], [-1, 0, 1], format='csc')
+
+    ug[k, :] = scipy.sparse.linalg.spsolve(A, b)
+
+  return ug

src/solvers/theoretical.py

@@ -0,0 +1,35 @@
+import functools
+import pickle
+
+import numpy as np
+
+from util import Grid
+
+with open('./roots128.pickle', 'rb') as file:
+  roots128 = pickle.load(file)
+
+
+def u_elem128(x, t, alpha_n):
+  x = np.float128(x)
+  t = np.float128(t)
+  alpha_n = np.float128(alpha_n)
+  return -4 * np.pi ** 2 * (alpha_n * np.cos(alpha_n * x / 10) + 40 * np.sin(alpha_n * x / 10)) * (
+      alpha_n * np.sin(2 * alpha_n) - 40 * np.cos(2 * alpha_n) + 40) * np.exp(-3 * alpha_n ** 2 * t / 500) / (
+             5 * (alpha_n ** 2 - np.pi ** 2) * (alpha_n ** 2 * (4 * alpha_n + np.sin(4 * alpha_n)) - 80 * alpha_n * (
+             np.cos(4 * alpha_n) - 1) + 6400 * alpha_n - 1600 * np.sin(4 * alpha_n)))
+
+
+@functools.lru_cache(maxsize=1000000)
+def u128(x, t, root_n=len(roots128)):
+  return sum(u_elem128(x, t, roots128[:root_n]))
+
+
+#@cachier(stale_after=datetime.timedelta(hours=1))
+def solve_theoretical(grid: Grid) -> np.ndarray:
+  ug = np.zeros((grid.K + 1, grid.I + 1), dtype=float)
+
+  for k, t in enumerate(grid.ts):
+    for i, x in enumerate(grid.xs):
+      ug[k, i] = u128(x, t, 1000)
+
+  return ug

src/util.py

@@ -0,0 +1,45 @@
+from dataclasses import dataclass, field
+import numpy as np
+
+import consts
+
+
+@dataclass
+class Grid:
+  I: int
+  K: int
+
+  xs: np.ndarray = field(repr=False, compare=False)
+  ts: np.ndarray = field(repr=False, compare=False)
+
+  h_x: float = field(repr=False, compare=False)
+  h_t: float = field(repr=False, compare=False)
+
+  gamma: float = field(repr=False, compare=False)
+
+  def __init__(self, I: int, K: int):
+    self.I = I
+    self.K = K
+
+    self.xs = np.linspace(0, consts.l, I + 1)
+    self.ts = np.linspace(0, consts.T, K + 1)
+
+    self.h_x = consts.l / I
+    self.h_t = consts.T / K
+
+    self.gamma = consts.D * self.h_t / (self.h_x ** 2)
+
+  def __hash__(self):
+    return hash((self.I, self.K))
+
+
+def mse(ug1: np.ndarray, ug2: np.ndarray) -> float:
+  return np.mean((ug1 - ug2) ** 2)
+
+
+def mae(ug1: np.ndarray, ug2: np.ndarray) -> float:
+  return np.max(np.abs(ug1 - ug2))
+
+
+def mae_last(ug1: np.ndarray, ug2: np.ndarray) -> float:
+  return np.max(np.abs(ug1[-1, :] - ug2[-1, :]))