1   %{
  2   Implementation of Adaptive Gauss-Kronrod Quadrature (GK15,7) 
  3   Octave 9.1.0 
  4   John Bryan 2025
  5   %}
  6 
  7 
  8   close all;
  9   clear all;
 10   echo off;
 11   clc;
 12   format long;
 13 
 14 
 15   function u=gk(f,a,b)
 16     % K15  
 17 
 18     x=[-0.9914553711208126,
 19        -0.9491079123427585,
 20        -0.8648644233597691,
 21        -0.7415311855993944,
 22        -0.5860872354676911,
 23        -0.4058451513773972,
 24        -0.2077849550078985,
 25         0.0000000000000000,
 26         0.2077849550078985,
 27         0.4058451513773972,
 28         0.5860872354676911,
 29         0.7415311855993944,
 30         0.8648644233597691,
 31         0.9491079123427585,
 32         0.9914553711208126];
 33 
 34      w=[0.0229353220105292,
 35         0.0630920926299786,
 36         0.1047900103222501,
 37         0.1406532597155259,
 38         0.1690047266392679,
 39         0.1903505780647854,
 40         0.2044329400752988,
 41         0.2094821410847278,
 42         0.2044329400752988,
 43         0.1903505780647854,
 44         0.1690047266392679,
 45         0.1406532597155259,
 46         0.1047900103222501,
 47         0.0630920926299786,
 48         0.0229353220105292];
 49 
 50     c=(b-a)/2;
 51     d=(b+a)/2;
 52     sum0=0;
 53     for k=[1:1:15]
 54         sum0=sum0+w(k)*f(c*x(k)+d);
 55     endfor
 56     u=c*sum0;
 57     return;
 58   endfunction
 59 
 60 
 61   function u1=gl(f,a,b)
 62           % G(7)
 63 
 64     x1=[-0.94910791234275852452618968404785126,
 65          -0.7415311855993944398638647732807884,
 66          -0.4058451513773971669066064120769615,
 67           0.0000000000000000000000000000000000,
 68           0.4058451513773971669066064120769615,
 69           0.7415311855993944398638647732807884,
 70           0.9491079123427585245261896840478513];
 71 
 72      w1=[0.1294849661688696932706114326790820,
 73          0.2797053914892766679014677714237796,
 74          0.3818300505051189449503697754889751,
 75          0.4179591836734693877551020408163265,
 76          0.3818300505051189449503697754889751,
 77          0.2797053914892766679014677714237796,
 78          0.1294849661688696932706114326790820];
 79 
 80     c=(b-a)/2;
 81     d=(b+a)/2;
 82     sum1=0;
 83     for k=[1:1:7]
 84         sum1=sum1+w1(k)*f(c*x1(k)+d);
 85     endfor
 86     u1=c*sum1;
 87     return;
 88   endfunction
 89 
 90 
 91   function [interval_result,error0] = adaptive_gk(f,a,b,abstol)
 92         # recursive adaptive gk 
 93         # f: integrand
 94         # a: lower limit
 95         # b: upper limit
 96         # abstol: error abstol
 97         s1 = gl(f,a,b);
 98         s2 = gk(f,a,b);
 99         error0 = abs(s1-s2);
100         if (error0 < abstol)
101            interval_result = s2;
102         else
103            [ltotal, error_ltotal] = adaptive_gk(f,a,(a+b)/2,abstol/2);
104            [rtotal, error_rtotal] = adaptive_gk(f,(a+b)/2,b,abstol/2);
105            interval_result= ltotal+rtotal;
106            error0 = error_ltotal+error_rtotal;
107         endif
108   endfunction
109 
110 
111   function test()
112       a=0;
113       b=50;
114       f=@(x) sin(x);
115       abstol =10^(-14);
116       [result,error0] = adaptive_gk(f,a,b,abstol);
117       exact_value=-cos(b)+cos(a);
118       difference=abs(result-exact_value);
119       fprintf("\nsin(x), 0:50 : Adaptive GK result %1.15e  error %1.15e\n\n",result,difference);
120       fprintf("eps = %1.15e\n\n",eps);
121       % the error is eps
122   endfunction
123 
124 
125 test()