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    # GK15
 13    # output GK15 quadrature result 
 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      # GL(7)
 57      # output GL quadrature result 
 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     # recursive adaptive gk 
 86     # f: integrand
 87     # a: lower limit
 88     # b: upper limit
 89     # abstol: error abstol
 90     # output interval result, error
 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()