Merge pull request #119 from eric2003/hexin

add eno example

Author: eric2003

example/1d-linear-convection/analytical_solution/01/theory.py

@@ -0,0 +1,54 @@
+import numpy as np
+
+def initial_condition(x):
+    """生成初始条件：在 [0.1, 0.3] 区间内为 1.0，其他位置为 0.0。
+    
+    参数:
+        x (np.ndarray): 空间坐标数组
+    
+    返回:
+        np.ndarray: 初始条件数组
+    """
+    return np.where((x >= 0.1) & (x <= 0.3), 1.0, 0.0)
+
+def analytical_solution(x, t, a, domain_length):
+    """计算理论解：初始波形以速度 a 平移，并应用周期边界条件。
+    
+    参数:
+        x (np.ndarray): 空间坐标数组
+        t (float): 时间
+        a (float): 波速
+        domain_length (float): 区域的周期长度 L
+    
+    返回:
+        np.ndarray: 时间 t 后的波形
+    """
+    # 计算平移后的坐标（无边界处理）
+    shifted_x = x - a * t
+    
+    # 应用周期边界条件，将坐标限制在 [0, domain_length) 范围内
+    shifted_x_periodic = (shifted_x + domain_length) % domain_length
+    
+    # 返回平移后的初始条件
+    return initial_condition(shifted_x_periodic)
+    
+    
+import matplotlib.pyplot as plt
+
+# 定义参数
+L = 1.0       # 区域长度
+a = 1.0       # 波速
+t = 0.4       # 时间
+x = np.linspace(0, L, 1000)  # 空间网格
+
+# 计算初始条件和理论解
+u0 = initial_condition(x)
+u_analytical = analytical_solution(x, t, a, L)
+
+# 可视化
+plt.plot(x, u0, label="Initial condition")
+plt.plot(x, u_analytical, label=f"Analytical (t={t})")
+plt.legend()
+plt.xlabel("x")
+plt.ylabel("u")
+plt.show()    

example/1d-linear-convection/analytical_solution/test_periodic_mapping/01/test_periodic_mapping.py

@@ -0,0 +1,18 @@
+import numpy as np
+
+def test_periodic_mapping():
+    L = 1.0
+    test_cases = [
+        (-0.2, 0.8),
+        (1.3, 0.3),
+        (0.5, 0.5),
+        (-1.1, 0.9),
+        (2.7, 0.7),
+    ]
+    
+    for x, expected in test_cases:
+        result = (x + L) % L
+        assert np.isclose(result, expected), f"Failed for {x}: {result} vs {expected}"
+    print("所有测试通过！")
+
+test_periodic_mapping()

example/1d-linear-convection/eno/eno_coef/01/crj.py

@@ -0,0 +1,46 @@
+import numpy as np
+
+def calculate_eno_crj(r, j, k):
+    result = 0
+
+    # 外层求和：m 从 j+1 到 k
+    for m in range(j + 1, k + 1):
+        numerator = 0
+
+        # 内层求和：l 从 0 到 k，且 l ≠ m
+        for l in range(0, k + 1):
+            if l == m:
+                continue  # 跳过 l = m 的情况
+
+            product = 1
+
+            # 内层连乘：q 从 0 到 k，且 q ≠ m, l
+            for q in range(0, k + 1):
+                if q == m or q == l:
+                    continue  # 跳过 q = m 或 q = l 的情况
+                product *= (r - q + 1)
+
+            numerator += product
+
+        denominator = 1
+
+        # 分母连乘：l 从 0 到 k，且 l ≠ m
+        for l in range(0, k + 1):
+            if l == m:
+                continue  # 跳过 l = m 的情况
+            denominator *= (m - l)
+
+        result += numerator / denominator
+
+    return result
+    
+for k in range(1,8):
+    print(f"=== k = {k} ===")
+    # 计算矩阵并存储到列表中
+    eno_mat = np.zeros((k+1,k))
+    for r in range(-1,k):
+        for j in range(k):
+            eno_mat[r+1, j] = calculate_eno_crj(r, j, k)
+        
+    print(f'{eno_mat=}')
+    print(f'{type(eno_mat)=}')