数值积分——牛顿-柯特斯公式

%%牛顿-柯特斯公式（此方法对于8阶以下是有效的，8阶以上误差将非常大）
%%interva为求积区间，Y随attribute变化（0或1）而对应不同选项（已知X对应的数值 或 表达式），n为步数
function NCF = Newton_Cotes_Formula(interval,Y,attribute,n)
a = interval(1);b = interval(2);
h = (b-a)/n;
for i = 1:n+1
    X(i) = a+h*(i-1);
end
if attribute == 0
    [m1,n1] = size(Y);
    MAX = max([m1,n1]);
    if MAX < n+1
        p = ceil((n+1)/MAX);
        Y_charge(:,1) = reshape(Y(1:MAX-1),MAX-1,1);
        lambda = input(‘输入插值因子（介于0,1之间）：‘);
        for i = 1:MAX-1
            for k = 2:p
                Y_charge(i,k) = Y(i)*(k-1)*(lambda/p)+Y(i+1)*(1-(k-1)*(lambda/p));
            end
        end
        Y_charge0 = reshape(Y_charge‘,1,(n1-1)*p);
        r = rand(1);
        Y_charge0(1,(n1-1)*p+1) = lambda*r*Y_charge0((n1-1)*p)+(1-r*lambda)*Y(MAX);
        Y = [Y_charge0,Y(MAX)];
        Y = Y(1:n+1);
    elseif MAX >n+1
        for i = 1:n+1
            X_charge(i) = floor(i*MAX/(n+1));
            Y_charge(i) = Y(X_charge(i));
        end
        Y = Y_charge;
    end
    sum = 0;
    for k = 1:n+1
        sum = sum+Cotes_coefficient(k-1,n)*Y(k);
    end
    NCF = vpa((b-a)*sum,6);
elseif attribute == 1
    a = interval(1);b = interval(2);
    h = (b-a)/n;
    for i = 1:n+1
        X(i) = a+h*(i-1);
    end
    F = subs(Y,X);
    sum = 0;
    for k = 1:n+1
        sum = sum+Cotes_coefficient(k-1,n)*F(k);
    end
    NCF = vpa((b-a)*sum,6);
end
%%柯特斯系数
    function CC = Cotes_coefficient(k,n)
        t = sym(‘t‘);
        mult = 1;
        for j = 0:1:n
            if j ~=k
                mult = mult*(t-j);
            end
        end
        C = int(mult,0,n);
        CC = (-1)^(n-k)/(n*factorial(k)*factorial(n-k))*C;
    end

    function F = factorial(n)
        if n == 0
            F = 1;
        else
            F = factorial(n-1)*n;
        end
    end

end

syms x;
Y = exp(x)*sin(x)+log(x+1);
interval=[0 pi];
attribute = 1;
n = 6;
Newton_Cotes_Formula(interval,Y,attribute,n)

vpa(int(Y,x,interval),6)

ans =
14.8156
ans =
14.8143