Project 5

文章发布时间:

最后更新时间:

练习一：Simpson 求积公式

实验发现，Simpson 方法求积小数点后只有一位有效数字，而采用复合Simpson 方法小数点后有五位有效数字。

代码： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
# Simpson 求积公式
from cmath import pi
import scipy as sc
import numpy as np

def f(x):
    M = 200
    sig = 0.1
    xbar = 1.7
    return M/(sig*(2*pi)**(1/2))*np.exp(-(x-xbar)**2/(2*sig**2))

def simpson(a,b):
    I = (b-a)/6*(f(a)+f(b)+4*f((a+b)/2))
    return I

a = 1.8
b = 1.9

I1 = sc.integrate.quad(f,a,b)
print(I1)

I2 = simpson(a,b)
print(I2)

I3 = 0
n = 10
h = (b-a)/n
for i in range(n):
    a1 = a + i*h
    b1 = a + (i+1)*h
    I3 = I3 + simpson(a1,b1)
print(I3)
输出结果： 
1
2
3
(27.181024396655484, 3.0176999118065023e-13) # 精确值
27.134402456529873	# Simpson 求积公式
27.18102028106211 # 复合Simpson 求积

练习二：Romberg 求积方法)

实验发现，取公差为0.000001，代码计算得到的近似值与“精确值”的小数点后七位都是相同的。

代码： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
# Romberg 求积方法
from cmath import pi
import scipy as sc
import numpy as np

def trapezoidal(a,b):  # 梯形求积
    I = (b-a)/2*(f(a)+f(b))
    return I

def romberg(f,a,b,eps):
    e = 1
    m = 2  # 初始阶
    while e > eps:
        h = b-a
        T = np.zeros((m,m))
        for k in range(m):
            n = np.power(2,k)
            h1 = h/n
            I = 0
            for s in range(n): #复合求积公式
                a1 = a + s*h1
                b1 = a + (s+1)*h1
                I = I + trapezoidal(a1,b1)
            T[k,0] = I
            for j in range(1,k+1):
                T[k,j]=1/(np.power(4,j)-1)*(np.power(4,j)*T[k,j-1]-T[k-1,j-1])
                #print(T) #测试
        e = abs(T[-1,-1]-T[-1,-2]) # 公差
        m += 1
        #print(T) #测试
    return T[-1,-1]

f = lambda x: np.log(x)
a,b = 1,2
eps = 0.000001
I1 = sc.integrate.quad(f,a,b)
print(I1[0])

I2 = romberg(f,a,b,eps)
print(I2)
输出结果： 
1
2
0.38629436111989063 #精确值
0.38629430908624807 #Romberg 算法

练习三：复合积分方法和 Gauss 积分方法的比较

复合方法求积

代码： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
# 复合积分方法和 Gauss 积分方法的比较
# 复合积积分方法
from cmath import pi
import scipy as sc
import numpy as np

def f(x):
    return np.sqrt(1+np.exp(x))

def f1(x):
    return 2*np.sqrt(1+np.exp(2*x+2))

def simpson(a,b):
    I = (b-a)/6*(f(a)+f(b)+4*f((a+b)/2))
    return I

def trapezoidal(a,b):  # 梯形求积
    I = (b-a)/2*(f(a)+f(b))
    return I


a = 0
b = 4

I1 = sc.integrate.quad(f,a,b)
print("numpy计算结果:",I1)

I3 = 0
I4 = 0
#n = 10
N = [2,4,8,16]
for n in N:
    h = (b-a)/n
    I3,I4=0,0
    for i in range(n):
        a1 = a + i*h
        b1 = a + (i+1)*h
        I3 = I3 + trapezoidal(a1,b1)
        I4 = I4 + simpson(a1,b1)
    print("复合梯形求积结果:",I3,"\n与精确解偏差:",I3-I1[0]) #误差以4倍速度下降 收敛阶为2
    print("复合Simpson求积结果",I4,"\n与精确解偏差:",I4-I1[0])#误差收敛阶为4
    print("\n")
输出结果： 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
numpy计算结果: (13.577302400789666, 1.5073833737655101e-13)
复合梯形求积结果: 14.663403727504809 
与精确解偏差: 1.0861013267151431
复合Simpson求积结果 13.581379562032346
与精确解偏差: 0.004077161242680205


复合梯形求积结果: 13.851885603400463
与精确解偏差: 0.2745832026107973
复合Simpson求积结果 13.57756442812153
与精确解偏差: 0.0002620273318640898


复合梯形求积结果: 13.646144721941265
与精确解偏差: 0.0688423211515996
复合Simpson求积结果 13.577318889241141
与精确解偏差: 1.6488451475282773e-05

复合梯形求积结果: 13.594525347416171
与精确解偏差: 0.017222946626505475
复合Simpson求积结果 13.577303433056535
与精确解偏差: 1.032266869316345e-06


使用Gauss节点结果(N= 3 ): 13.576984579216827
与精确解偏差: -0.00031782157283899437
使用Gauss节点结果(N= 5 ): 13.57730258966976
与精确解偏差: 1.888800937877022e-07
使用Gauss节点结果(N= 7 ): 13.577302399995329
与精确解偏差: -7.943370405882888e-10
PS C:\Users\bc\OneDrive\桌面\数值分析> & C:/Users/bc/AppData/Local/Programs/Python/Python310/python.exe c:/Users/bc/OneDrive/桌面/数值分析/数 
值逼近/11.09.01.py
numpy计算结果: (13.577302400789666, 1.5073833737655101e-13)
复合梯形求积结果: 14.663403727504809 
与精确解偏差: 1.0861013267151431
复合Simpson求积结果 13.581379562032346
与精确解偏差: 0.004077161242680205


复合梯形求积结果: 13.851885603400463
与精确解偏差: 0.2745832026107973
复合Simpson求积结果 13.57756442812153
与精确解偏差: 0.0002620273318640898


复合梯形求积结果: 13.646144721941265
与精确解偏差: 0.0688423211515996
复合Simpson求积结果 13.577318889241141
与精确解偏差: 1.6488451475282773e-05


复合梯形求积结果: 13.594525347416171
与精确解偏差: 0.017222946626505475
复合Simpson求积结果 13.577303433056535
与精确解偏差: 1.032266869316345e-06
## Gauss 积分公式求积 代码： 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
# 复合积分方法和 Gauss 积分方法的比较
# 复合积积分方法
from cmath import pi
import scipy as sc
import numpy as np

def f(x):
    return np.sqrt(1+np.exp(x))

def f1(x):
    return 2*np.sqrt(1+np.exp(2*x+2))


a = 0
b = 4

I1 = sc.integrate.quad(f,a,b)
I3 = 0
I4 = 0

# Gauss 积分公式
N3=np.array([[0.8888888888888888,0.0000000000000000],
[0.5555555555555556,-0.7745966692414834],
[0.5555555555555556,0.7745966692414834]])

N5=np.array([[0.5688888888888889, 0.0000000000000000]
,[0.4786286704993665, -0.5384693101056831]
,[0.4786286704993665, 0.5384693101056831]
,[0.2369268850561891, -0.9061798459386640]
,[0.2369268850561891, 0.9061798459386640]])

N7=np.array([[0.4179591836734694, 0.0000000000000000],
[0.3818300505051189, -0.4058451513773972],
[0.3818300505051189, 0.4058451513773972],
[0.2797053914892766, -0.7415311855993945],
[0.2797053914892766, 0.7415311855993945],
[0.1294849661688697, -0.9491079123427585],
[0.1294849661688697, 0.9491079123427585]])

for N in [N3,N5,N7]:
    I5 = 0
    for i in N:
        I5 += 1*i[0]*f1(i[1])
    print("使用Gauss节点结果(N=",len(N),"):",I5)
    print("与精确解偏差:",I5-I1[0])
输出结果： 
1
2
3
4
5
6
使用Gauss节点结果(N= 3 ): 13.576984579216827
与精确解偏差: -0.00031782157283899437
使用Gauss节点结果(N= 5 ): 13.57730258966976
与精确解偏差: 1.888800937877022e-07
使用Gauss节点结果(N= 7 ): 13.577302399995329
与精确解偏差: -7.943370405882888e-10

≡