function [ integration ] = CompoInt( func, left, right, step, mode )
%   分段积分牛顿-柯茨公式；
%   输入5个参数：被积函数(FUNCTIONHANDLE)\'func\'，积分上下界\'left\'/\'right\'，步长\'step\'，
%   模式mode共三种(\'midpoint\',\'trapezoid\',\'simpson\')；
%   返回积分值；
switch mode
    % 仅仅是公式不同
    case \'midpoint\'
        node = left; integration = 0;
        while (node + step <= right)    % 按照给定步长开始分段积分
            pieceInt = step*(func(node + step/2));  % 使用中点积分公式
            integration = integration + pieceInt; node = node + step;
        end
        if (node < right)   % 补齐最后一段积分
            pieceInt = (right - node)*(func((node + right)/2));
            integration = integration + pieceInt;
        end
    case \'trapezoid\'
        node = left; integration = 0;
        while (node + step <= right)
            pieceInt = step*(func(node) + func(node + step))/2; % 使用梯形公式
            integration = integration + pieceInt; node = node + step;
        end
        if (node < right)
            pieceInt = (right - node)*(func(node) + func(right))/2;
            integration = integration + pieceInt;
        end
    case \'simpson\'
        node = left; integration = 0;
        while (node + step <= right)
            pieceInt = step*(func(node) + 4*func(node + step/2) + func(node + step))/6; % 使用辛普森公式
            integration = integration + pieceInt; node = node + step;
        end
        if (node < right)
            pieceInt = (right - node)*(func(node) + 4*func((node + right)/2) + func(right))/6;
            integration = integration + pieceInt;
        end
end

function [ integration ] = AdaptInt( func, left, right, prec, length )
%   自适应分段积分函数AdaptInt；
%   输入五个参数：被积函数（句柄）func，积分上下限right，left，要求精度prec，积分总长length；
%   输出一个参数：积分值integration；
trapeInt = (right - left)*(func(left) + func(right))/2;
midInt = (right - left)*func((left + right)/2);
err = (trapeInt - midInt)/3;    % 由中点公式和梯形公式差估算误差
if (abs(err) < prec && (right - left) < length/5)   % 如果误差符合要求，则使用辛普森公式计算较精确结果
    integration = (right - left)*(func(left) + 4*func((left + right)/2) + func(right))/6;
else    % 否则，二分该段，分别进行自适应分段积分
    integration = AdaptInt(func, left, (left + right)/2, prec/2, length) + AdaptInt(func, (left + right)/2, right, prec/2, length);
end
end