數值積分：基於牛頓-柯茨公式的定步長和自適應積分方法 [MATLAB]

#先上代碼后補筆記#

#可以直接復制粘貼使用的MATLAB函數！#

1. 定步長牛頓-柯茨積分公式

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

　　

2. 自適應分段積分方法

通過函數內調用自身的方法使得積分函數顯得簡潔明快。因為需要調用自身計算原積分的一段，必須傳入參數length標識原積分上下限總長，通過left, right和length三個參數才能夠得到某一段的要求精度。

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