Переписать трапецеидальное правило Симпсона в Matlab - PullRequest
Новая
       35

Переписать трапецеидальное правило Симпсона в Matlab

0 голосов
/ 01 мая 2018

Я пытаюсь написать программу Matlab для вычисления интеграла с помощью правила трапеции и Симпсонов. Программа для трапеции выглядит следующим образом: 

    function [int, flag, stats] = trapComp(f, a, b, tol, hMin)
        % Initialise variables
        h = b - a;
        n = 1;
        int = h / 2 * (f(a) + f(b));
        flag = 1;

        if nargout == 3
            stats = struct('totalErEst', [], 'totalNrIntervals', [], 'nodesList', []);
        end

        while h > hMin
            h = h / 2;
            n = 2 * n;
            if h < eps % Check if h is not "zero"
                break;
            end

            % Update the integral with the new nodes
            intNew = int / 2;
            for j = 1 : 2 : n
                intNew = intNew + h * f(a + j * h);
            end

            % Estimate the error
            errorEst = 1 / 3 * (int - intNew);
            int = intNew;
            if nargout == 3 % Update stats
                stats.totalErEst = [stats.totalErEst; abs(errorEst)];
                stats.totalNrIntervals = [stats.totalNrIntervals; n / 2];
            end

            if abs(errorEst) < tol
                flag = 0;
                break
            end
        end
    end

Теперь правила Симпсонов, я не могу обойтись. Я знаю, что это очень похоже, но я не могу понять это.

Это мой код Симпсона: 

function [int, flag, stats] = simpComp(f, a, b, tol, hMin)
    % Initialise variables
    h = b - a;
    n = 1;
    int = h / 3 * (f(a) + 4 * f((a+b)/2) + f(b));
    flag = 1;

    if nargout == 3
        stats = struct('totalErEst', [], 'totalNrIntervals', [], 'nodesList', []);
    end

    while h > hMin
        h = h / 2;
        n = 2 * n;
        if h < eps % Check if h is not "zero"
            break;
        end

        % Update the integral with the new nodes
        intNew = int / 2;
        for j = 1 : 2 : n
            intNew = intNew + h * f(a + j * h);
        end

        % Estimate the error
        errorEst = 1 / 3 * (int - intNew);
        int = intNew;
        if nargout == 3 % Update stats
            stats.totalErEst = [stats.totalErEst; abs(errorEst)];
            stats.totalNrIntervals = [stats.totalNrIntervals; n / 2];
        end

        if abs(errorEst) < tol
            flag = 0;
            break
        end
    end
end

Использование этого, однако, дает ответ для интеграла с большей ошибкой, чем трапецеидальный, который, я считаю, не должен.

Любая помощь будет оценена

...