1. 定步长牛顿-柯茨积分公式
function [ integration ] = CompoInt( func, left, right, step, mode )
% 分段积分牛顿-柯茨公式;
% 输入5个参数:被积函数(FUNCTIONHANDLE)'func',积分上下界'left'/'right',步长'step',
% 模式mode共三种('midpoint','trapezoid','simpson');
% 返回积分值;
switch mode
% 仅仅是公式不同
case 'midpoint'
node = left; integration = 0;
while (node + step <= right) % 按照给定步长开始分段积分
pieceInt = step*(func(node + step/2)); % 使用中点积分公式
integration = integration + pieceInt; node = node + step;
end
if (node < right) % 补齐最后一段积分
pieceInt = (right - node)*(func((node + right)/2));
integration = integration + pieceInt;
end
case 'trapezoid'
node = left; integration = 0;
while (node + step <= right)
pieceInt = step*(func(node) + func(node + step))/2; % 使用梯形公式
integration = integration + pieceInt; node = node + step;
end
if (node < right)
pieceInt = (right - node)*(func(node) + func(right))/2;
integration = integration + pieceInt;
end
case 'simpson'
node = left; integration = 0;
while (node + step <= right)
pieceInt = step*(func(node) + 4*func(node + step/2) + func(node + step))/6; % 使用辛普森公式
integration = integration + pieceInt; node = node + step;
end
if (node < right)
pieceInt = (right - node)*(func(node) + 4*func((node + right)/2) + func(right))/6;
integration = integration + pieceInt;
end
end
2. 自适应分段积分方法
通过函数内调用自身的方法使得积分函数显得简洁明快。因为需要调用自身计算原积分的一段,必须传入参数length标识原积分上下限总长,通过left, right和length三个参数才能够得到某一段的要求精度。
function [ integration ] = AdaptInt( func, left, right, prec, length )
% 自适应分段积分函数AdaptInt;
% 输入五个参数:被积函数(句柄)func,积分上下限right,left,要求精度prec,积分总长length;
% 输出一个参数:积分值integration;
trapeInt = (right - left)*(func(left) + func(right))/2;
midInt = (right - left)*func((left + right)/2);
err = (trapeInt - midInt)/3; % 由中点公式和梯形公式差估算误差
if (abs(err) < prec && (right - left) < length/5) % 如果误差符合要求,则使用辛普森公式计算较精确结果
integration = (right - left)*(func(left) + 4*func((left + right)/2) + func(right))/6;
else % 否则,二分该段,分别进行自适应分段积分
integration = AdaptInt(func, left, (left + right)/2, prec/2, length) + AdaptInt(func, (left + right)/2, right, prec/2, length);
end
end