码迷,mamicode.com
首页 > 其他好文 > 详细

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

时间:2020-01-30 17:27:31      阅读:145      评论:0      收藏:0      [点我收藏+]

标签:and   UNC   att   结果   不同   log   pre   end   elseif   

  此段代码是牛顿- 柯特斯数值积分法,代码如下:

  1.代码

%%牛顿-柯特斯公式(此方法对于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

  2.例子

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)

  3.结果

ans =
14.8156
ans =
14.8143

  第一项结果为牛顿-柯特斯数值积分计算结果,第二项为真实值

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

标签:and   UNC   att   结果   不同   log   pre   end   elseif   

原文地址:https://www.cnblogs.com/guliangt/p/12243008.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!