function?[sfcount]=acsimpson（fabtol）
%此方法為自適應(yīng)simpson方法求解定積分
%f為函數(shù)，a，b分別為積分區(qū)間，tol為精度
%q為積分值，fcount為計算函數(shù)值次數(shù)

%一下為準備輸入數(shù)據(jù)
f?=?fcnchk（f‘vectorized‘）;???????%定義內(nèi)聯(lián)函數(shù)
h=b-a;?????????????%步長
qu=[a?a+h/2?b];????%三個節(jié)點
qu（qu==0）=2^-1024;?%避免計算函數(shù)值出現(xiàn)NAN的情況
y=feval（fqu）;?????%計算函數(shù)節(jié)點處的值
s=y*[1?4?1]‘*h/6;??%Simpson積分
fcount=0;??????????%定義初始計數(shù)器
tol=tol*15;????????%調(diào)整精度

[sfcount]=subsim（fsquyfcounttol）;?%遞歸函數(shù)
fcount=fcount+3;???%計數(shù)器加上準備輸入數(shù)據(jù)的三次計算
end

%以下為遞歸函數(shù)
function?[sfcount]=subsim（fsquyfcounttol）
%以下為對qu二分，分別對兩個小區(qū)間求出Simpson積分
h=（qu（end）-qu（1））/2;
quadd=[qu（1）+h/2?qu（2）+h/2];
quadd（quadd==0）=2^-1024;?%避免計算函數(shù)值出現(xiàn)NAN的情況
yadd=feval（fquadd）;
subs=[y（1）?yadd（1）?