1 # Adaptive Simpson
 2 # Octave 9.1.0
 3 # John Bryan 2024
 4 clear all;
 5 close all;
 6 clc;
 7 format long;
 8 echo off;
 9 
10 
11 function [interval_result,error0] = adaptive_simpson(f,a,b,abstol)
12     # recursive adaptive simpson
13     # f: integrand
14     # a: lower limit
15     # b: upper limit
16     # abstol: error abstol
17     h1=(b-a)/2;
18     h2 = h1/2;
19     s1 = h1/3*(f(a)+4*f(a+h1)+f(b));
20     s2 = h2/3*(f(a)+4*f(a+h2)+2*f(a+2*h2)+4*f(a+3*h2)+f(b));
21     error0 = abs(s1-s2)/15;
22     if (error0 < abstol)
23        interval_result = s2;
24     else
25        [ltotal, error_ltotal] = adaptive_simpson(f,a,(a+b)/2,abstol/2);
26        [rtotal, error_rtotal] = adaptive_simpson(f,(a+b)/2,b,abstol/2);
27        interval_result= ltotal+rtotal;
28        error0 = error_ltotal+error_rtotal;
29     endif
30 endfunction
31 
32 
33 function test()
34     f=@(x) sin(x);
35     a=0;
36     b=10;
37     abstol = eps;
38     exact_value=-cos(b)+cos(a);
39     [result,error0] = adaptive_simpson(f,a,b,abstol);
40     s4difference=abs(result-exact_value);
41     fprintf("\nsin(x), (%d,%d) : Adaptive Simpson result %1.15e  error %1.15e\n\n",a,b,result,s4difference);
42 endfunction
43 
44 
45 test()