1
2 '''
3 Adaptive Gauss Kronrod (GK15,7) Quadrature
4 Python 3.10.12
5 John Bryan, 2024
6 '''
7
8 import numpy as np
9
10
11 def gk(f,a,b):
12
13
14 xo=[-0.9914553711208126,
15 -0.9491079123427585,
16 -0.8648644233597691,
17 -0.7415311855993944,
18 -0.5860872354676911,
19 -0.4058451513773972,
20 -0.2077849550078985,
21 -0.0000000000000000,
22 0.2077849550078985,
23 0.4058451513773972,
24 0.5860872354676911,
25 0.7415311855993944,
26 0.8648644233597691,
27 0.9491079123427585,
28 0.9914553711208126];
29
30 wo=[0.0229353220105292,
31 0.0630920926299786,
32 0.1047900103222501,
33 0.1406532597155259,
34 0.1690047266392679,
35 0.1903505780647854,
36 0.2044329400752988,
37 0.2094821410847278,
38 0.2044329400752988,
39 0.1903505780647854,
40 0.1690047266392679,
41 0.1406532597155259,
42 0.1047900103222501,
43 0.0630920926299786,
44 0.0229353220105292];
45
46 c=(b-a)/2;
47 d=(b+a)/2;
48 sum0=0.;
49 for k in range(0,15):
50 sum0=sum0+wo[k]*f(c*xo[k]+d);
51 u=c*sum0;
52 return u
53
54
55 def gl(f,a,b):
56
57
58
59 x1=[-0.94910791234275852452618968404785126,
60 -0.7415311855993944398638647732807884,
61 -0.4058451513773971669066064120769615,
62 0.0000000000000000000000000000000000,
63 0.4058451513773971669066064120769615,
64 0.7415311855993944398638647732807884,
65 0.9491079123427585245261896840478513];
66
67 w1=[0.1294849661688696932706114326790820,
68 0.2797053914892766679014677714237796,
69 0.3818300505051189449503697754889751,
70 0.4179591836734693877551020408163265,
71 0.3818300505051189449503697754889751,
72 0.2797053914892766679014677714237796,
73 0.1294849661688696932706114326790820];
74
75 c=(b-a)/2;
76 d=(b+a)/2;
77 sum1=0.;
78 for k in range(0,7):
79 sum1=sum1+w1[k]*f(c*x1[k]+d);
80 u1=c*sum1;
81 return u1
82
83
84 def adaptive_gk(f,a,b,abstol):
85
86
87
88
89
90
91 s1 = gl(f,a,b)
92 s2 = gk(f,a,b)
93 error0 = abs(s1-s2)
94 if (error0 < abstol):
95 interval_result = s2
96 else:
97 [rtotal, error_rtotal] = adaptive_gk(f,(a+b)/2,b,abstol/2)
98 [ltotal, error_ltotal] = adaptive_gk(f,a,(a+b)/2,abstol/2)
99 interval_result= ltotal+rtotal
100 error0 = error_ltotal+error_rtotal
101 return interval_result,error0
102
103
104 def test():
105 a=0.
106 b=50.
107 f=np.sin
108 exact =-np.cos(50)+np.cos(0)
109 abstol=1e-14
110 [result,error0] = adaptive_gk(f,a,b,abstol);
111 difference=abs(result-exact)
112 print("The quadrature result is ", result)
113 print("The quadrature difference from exact is ", difference)
114
115
116 test()