复化梯形、复化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;endend

复化梯形求积

    计算公式如下：

   

    实现代码：

function [result] = ComplexTrap( x_LowBound, x_Up_Bound,n)% Inputs:%        x_LowBound:积分区间下界%        x_UpBound :积分区间上界%        n         ：等分数量% Outputs:%        result    : 复化梯形积分结果% 获取步长hstep_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 % 函数结束