forked from chfenger/goNum
-
Notifications
You must be signed in to change notification settings - Fork 0
/
IntegralGaussLagendre.go
130 lines (123 loc) · 3.76 KB
/
IntegralGaussLagendre.go
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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
// IntegralGaussLagendre
/*
------------------------------------------------------
作者 : Black Ghost
日期 : 2018-12-12
版本 : 0.0.0
------------------------------------------------------
不超过8次的Gauss-Lagendre求积分公式
理论：
对于积分
b
|f(x)dx
a
使用n+1次Lagendre多项式的零点作为高斯点，可获得代数
精度为2n+1的高斯型求积公式
其中区间[a, b]需要预先转换为区间[-1, 1]
参考 李信真, 车刚明, 欧阳洁, 等. 计算方法. 西北工业大学
出版社, 2000, pp 162-164.
------------------------------------------------------
输入 :
fun 被积分函数
a, b 积分区间
n 求积分公式次数
输出 :
sol 解
err 解出标志：false-未解出或达到步数上限；
true-全部解出
------------------------------------------------------
*/
package goNum
// IntegralGaussLagendre 不超过8次的Gauss-Lagendre求积分公式
func IntegralGaussLagendre(fun func(float64) float64, a, b float64, n int) (float64, bool) {
/*
不超过8次的Gauss-Lagendre求积分公式
输入 :
fun 被积分函数
a, b 积分区间
n 求积分公式次数
输出 :
sol 解
err 解出标志：false-未解出或达到步数上限；
true-全部解出
*/
//判断n范围
if (n < 1) || (n > 8) {
panic("Error in goNum.IntegralGaussLagendre: n is a not correct input")
}
xi := [][]float64{
{0.0},
{-0.5773502692, 0.5773502692},
{-0.7745966692, 0.0, 0.7745966692},
{-0.8611363116, -0.3399810436, 0.3399810436, 0.8611363116},
{-0.9061798459, -0.5384693101, 0.0, 0.5384693101, 0.9061798459},
{-0.9324695142, -0.6612093865, -0.2386191861, 0.2386191861, 0.6612093865, 0.9324695142},
{-0.9491079123, -0.7415311856, -0.4058451514, 0.0, 0.4058451514, 0.7415311856, 0.9491079123},
{-0.9602898566, -0.7966664774, -0.5255324099, -0.1834346425, 0.1834346425, 0.5255324099, 0.7966664774, 0.9602898566},
}
Ai := [][]float64{
{2.0},
{1.0, 1.0},
{0.555555555555556, 0.888888888888889, 0.555555555555556},
{0.3478548451, 0.6521451549, 0.6521451549, 0.3478548451},
{0.2369268851, 0.4786286705, 0.568888889, 0.4786286705, 0.2369268851},
{0.1713244924, 0.3607615730, 0.4679139346, 0.4679139346, 0.3607615730, 0.1713244924},
{0.1294849662, 0.2797053915, 0.3818300505, 0.4179591837, 0.3818300505, 0.2797053915, 0.1294849662},
{0.1012285363, 0.2223810345, 0.3137066459, 0.3626837834, 0.3626837834, 0.3137066459, 0.2223810345, 0.1012285363},
}
//区间转换
c := (b - a) / 2.0
d := (a + b) / 2.0
switch n {
case 1:
sol := 0.0
for i := 0; i < len(xi[0]); i++ {
sol += Ai[0][i] * fun(d+c*xi[0][i])
}
return c * sol, true
case 2:
sol := 0.0
for i := 0; i < len(xi[1]); i++ {
sol += Ai[1][i] * fun(d+c*xi[1][i])
}
return c * sol, true
case 3:
sol := 0.0
for i := 0; i < len(xi[2]); i++ {
sol += Ai[2][i] * fun(d+c*xi[2][i])
}
return c * sol, true
case 4:
sol := 0.0
for i := 0; i < len(xi[3]); i++ {
sol += Ai[3][i] * fun(d+c*xi[3][i])
}
return c * sol, true
case 5:
sol := 0.0
for i := 0; i < len(xi[4]); i++ {
sol += Ai[4][i] * fun(d+c*xi[4][i])
}
return c * sol, true
case 6:
sol := 0.0
for i := 0; i < len(xi[5]); i++ {
sol += Ai[5][i] * fun(d+c*xi[5][i])
}
return c * sol, true
case 7:
sol := 0.0
for i := 0; i < len(xi[6]); i++ {
sol += Ai[6][i] * fun(d+c*xi[6][i])
}
return c * sol, true
case 8:
sol := 0.0
for i := 0; i < len(xi[7]); i++ {
sol += Ai[7][i] * fun(d+c*xi[7][i])
}
return c * sol, true
default:
return 0.0, false
}
}