91 lines
1.9 KiB
Python
91 lines
1.9 KiB
Python
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()
|