作业2 完成

numerical_analysis/2/README.md

@@ -0,0 +1,17 @@
+# 结论
+
+图例
+
+    0 - 标准函数
+    1 - 复化梯形公式
+    2 -  复化Simpson 1/3
+    3 -  复化Simpson 3/8
+
+随积分区间变化图
+
+![image-20210117003720172](https://public.veypi.com/img/screenshot/20210117003720.png)
+
+
+随N变化图
+
+![image-20210117004402891](https://public.veypi.com/img/screenshot/20210117004402.png)

numerical_analysis/2/main.py

@@ -0,0 +1,90 @@
+import math
+from scipy import integrate
+import matplotlib.pyplot as plt
+
+
+def erf_inner(x):
+    return math.e ** (- x ** 2)
+
+
+def solve0(a, b, N):
+    """
+    标准函数求积分
+    """
+    return integrate.quad(erf_inner, a, b)[0]
+
+
+def solve1(a, b, N):
+    """
+    复化梯形积分公式
+    """
+    h = (b - a) / N
+    s = erf_inner(a)
+    for i in range(1, N):
+        s += 2 * erf_inner(a + h * i)
+    s += erf_inner(b)
+    return s * h / 2
+
+
+def solve2(a, b, N):
+    """
+    复化Simpson 1/3
+    """
+    h = (b - a) / N
+    s = 0
+    s = erf_inner(a)
+    for i in range(1, N):
+        s += 2 * erf_inner(a + h * i)
+    for i in range(0, N):
+        s += 4 * erf_inner(a + h * (i + 0.5))
+    s += erf_inner(b)
+    return s * h / 6
+
+
+def solve3(a, b, N):
+    """
+    复化Simpson 3/8
+    """
+    h = (b - a) / N
+    s = 0
+    s = 7 * erf_inner(a)
+    for i in range(1, N):
+        s += 32 * erf_inner(a + h * (i - 0.75)) + 12 * erf_inner(a + h * (i - 0.5)) + 32 * erf_inner(
+            a + h * (i - 0.25)) + 14 * erf_inner(a + h * i)
+    s += 7 * erf_inner(b)
+    return s * h / 90
+
+
+def show():
+    x = [i for i in range(60)]
+    y = [[0 for j in range(60)] for _ in range(4)]
+    solves = [solve0, solve1, solve2, solve3]
+    for s in range(4):
+        for i in x:
+            y[s][i] = solves[s](0, i, 10)
+        plt.plot(x, y[s], label=str(s))
+    plt.legend(loc="lower right")
+    plt.show()
+
+
+def showN():
+    x = [i for i in range(1, 30)]
+    y = [[0 for j in range(30)] for _ in range(4)]
+    solves = [solve0, solve1, solve2, solve3]
+    for s in range(4):
+        for i in x:
+            y[s][i] = solves[s](0, 10, i)
+        plt.plot(x, y[s][1:], label=str(s))
+    plt.legend(loc="lower right")
+    plt.show()
+
+
+if __name__ == '__main__':
+    # a = 0
+    # b = 20
+    # n = 90
+    # print(solve1(a, b, n))
+    # print(solve2(a, b, n))
+    # print(solve3(a, b, n))
+    # print(solve0(a, b, n))
+    showN()