复化梯形、复化simpson以及复化Cotes的Matlab实现

复化求积法

    将区间[a,b]分成n等分，节点为xk = a + kh ( k = 0,1,2,....,n ), 步长 h = (b-a)/n,在每个子区间[xk,xk+1]上先用低阶Newton-Cotes公式，求得积分近似值Ik，再对k = 0,1,...,n-1对Ik求和作为准确积分I的近似值。

被积分函数

    此处以f(x)=sinx/x作为被积分函数进行说明，由于f(x)原函数很难找到，因此可以采用数值积分的方式进行计算。

    f(x)实现的代码为：

%计算x对应的函数值f(x)
function [y_x] = CalcuFunctionValue(x)
% inputs:
%        x:待求值
% outputs:
%        y_x:x对应的函数值

% 根据极限计算，当x→0时，y_x→1.
if x == 0
    y_x = 1;
else
    y_x = sin(x)/x;
end
end

复化梯形求积

    计算公式如下：

   

    实现代码：

function [result] = ComplexTrap( x_LowBound, x_Up_Bound,n)
% Inputs:
%        x_LowBound:积分区间下界
%        x_UpBound :积分区间上界
%        n         ：等分数量
% Outputs:
%        result    : 复化梯形积分结果

% 获取步长h
step_length = (x_Up_Bound - LowBound)/n;
%累积计算
result = 0;
for i = 1:n-1
    result = result + CalcuFunctionValue(x_LowBound+step_length*(i-1))+CalcuFunctionValue(x_LowBound+step_length*i);
end % 循环结束
result = result * step_length / 2;
end % 函数结束

复化Simpson求积

    计算公式如下：

   

    实现代码：

function [result] = ComplexSimpson(x_LowBound, x_Up_Bound,n)
% simpson求积公式
% Inputs:
%        x_LowBound:积分区间下界
%        x_UpBound :积分区间上界
%        n         ：等分数量,需要为2n等分，即节点个数必须满足2n+1
% Outputs:
%        result    : 复化Simpson积分结果

% 判断积分区间个数是否是2的倍数，满足则进行计算，否则打印提示
if mod(n,2) == 0
    % 获取步长h
    step_length = (x_Up_Bound - LowBound)/n;
    %累积计算
    result = 0;
    for i = 1:2:n-1
        result = result + CalcuFunctionValue(x_LowBound+step_length*(i-1))...
                 +4*CalcuFunctionValue(x_LowBound+step_length*i)...
                 +CalcuFunctionValue(x_LowBound+step_length*(i+1));
    end % 循环结束
    result = result*step_length/6;
else
    print('等分区间数错误！');
end % if判断结束
end % 函数结束

复化Cotes积分

    计算公式如下：

    

    代码实现：

%%
function [result] = ComplexCotes(x_LowBound, x_Up_Bound,n)
% cotes求积公式
% Inputs:
%        x_LowBound:积分区间下界
%        x_UpBound :积分区间上界
%        n         ：等分数量,需要为4n等分，即节点个数必须满足4n+1
% Outputs:
%        result    : 复化cotes积分结果

% 判断积分区间个数是否是2的倍数，满足则进行计算，否则打印提示
if mod(n,4) == 0
    % 获取步长h
    step_length = (x_Up_Bound - LowBound)/n;
    %累积计算
    result = 0;
    for i = 1:4:n-3
        result = result + 7*CalcuFunctionValue(x_LowBound+step_length*(i-1))...
                 +32*CalcuFunctionValue(x_LowBound+step_length*i)...
                 +12*CalcuFunctionValue(x_LowBound+step_length*(i+1))...
                 +32*CalcuFunctionValue(x_LowBound+step_length*(i+2))...
                 +7*CalcuFunctionValue(x_LowBound+step_length*(i+3));
    end % 循环结束
    result = result*step_length/90;
else
    print('等分区间数错误！');
end % if判断结束
end % 函数结束