1 # Fejer First Rule Adaptive Quadrature
 2 # Octave 9.1.0
 3 # John Bryan, 2025
 4 
 5 close all;
 6 clear all;
 7 echo off;
 8 clc;
 9 format long;
10 #syms x;
11 
12 
13 function u1=f0(f,a,b)
14    # x1: Chebyshev zeros: for k=0:1:(n-1) x1(k+1)=cos(((k+.5)*pi)/n); endfor
15    x1=[-0.98480775301220805936674302459,
16        -0.86602540378443864676372317075,
17        -0.64278760968653932632264340991,
18        -0.34202014332566873304409961468,
19         0.00000000000000000000000000000,
20         0.34202014332566873304409961468,
21         0.64278760968653932632264340991,
22         0.86602540378443864676372317075,
23         0.98480775301220805936674302459];
24    # w1: weights: w1(i)=integral from -1 to 1 of lagrange polynomial
25    #                                             Prod(x-x1(j))/Prod(x1(i)-x1(j)), j=[1:1:9], j!=i 
26    #                                                                           i=[1:1:9]
27    w1=[0.05273664990990678164016503879,
28        0.17918871252204585160012268514,
29        0.26403722254100439718081377617,
30        0.33084517516813642278027841712,
31        0.34638447971781303808608893430,
32        0.33084517516813642278027841712,
33        0.26403722254100439718081377617,
34        0.17918871252204585160012268514,
35        0.05273664990990678164016503879];
36    c=(b-a)/2;
37    d=(b+a)/2;
38    sum1=0;
39       for k=[1:1:9]
40          sum1=sum1+w1(k)*f(c*x1(k)+d);
41       endfor
42    u1=c*sum1;
43    return;
44 endfunction
45 
46 
47 function [interval_result,error0] = adaptive_fejer_1st_rule(f,a,b,abstol)
48    # f: integrand
49    # a: lower limit
50    # b: upper limit
51    # abstol: error abstol
52    s1 = f0(f,a,b);
53    s2 = f0(f,a,(a+b)/2);
54    s3 = f0(f,(a+b)/2,b);
55    s4 = s2+s3;
56    error0 = abs(s1-s4);
57    if (error0 < abstol)
58       interval_result = s1;
59    else
60       [ltotal, error_ltotal] = adaptive_fejer_1st_rule(f,a,(a+b)/2,abstol/2);
61       [rtotal, error_rtotal] = adaptive_fejer_1st_rule(f,(a+b)/2,b,abstol/2);
62       interval_result= ltotal+rtotal;
63       error0 = error_ltotal+error_rtotal;
64    endif
65 endfunction
66 
67 
68 function test()
69    n=9;
70    a=0;
71    b=50;
72    f=@(x) sin(x);
73    abstol=10^(-14);
74    [result,error0] = adaptive_fejer_1st_rule(f,a,b,abstol);
75    exact_value=-cos(b)+cos(a);
76    difference=abs(result-exact_value);
77    fprintf("\nsin(x), 0:50 : Adaptive Fejer First Rule Quadrature result %1.15e  error %1.15e\n\n",result,difference);
78    fprintf("eps = %1.15e\n\n",eps);
79    fprintf("the error is eps\n\n");
80 endfunction
81 
82 
83 test()