From 0393893d3dff54cd86e8a7274a23ee5df3618f70 Mon Sep 17 00:00:00 2001
From: veypi <veypi@qq.com>
Date: Sun, 17 Jan 2021 00:45:11 +0800
Subject: [PATCH] =?UTF-8?q?=E4=BD=9C=E4=B8=9A2=20=E5=AE=8C=E6=88=90?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 numerical_analysis/2/README.md | 17 +++++++
 numerical_analysis/2/main.py   | 90 ++++++++++++++++++++++++++++++++++
 2 files changed, 107 insertions(+)
 create mode 100644 numerical_analysis/2/main.py

diff --git a/numerical_analysis/2/README.md b/numerical_analysis/2/README.md
index e69de29..d959e58 100644
--- a/numerical_analysis/2/README.md
+++ b/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)
\ No newline at end of file
diff --git a/numerical_analysis/2/main.py b/numerical_analysis/2/main.py
new file mode 100644
index 0000000..7960b4d
--- /dev/null
+++ b/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()