我需要帮助找到使用梯形求和函数的积分。

程序应该采取连续梯形总和与n = 1, 2, 3, ...的子区间，直到有两个neighouring值n ，通过小于一个给定的容差有所不同。我想至少一个FOR一个内循环WHILE循环，我不希望使用trapz功能。该方案需要四个输入：

f ：用于的功能的功能句柄x 。
a ：实数。
b ：一个实数大于a 。
tolerance ：一个实数，它是积极的，非常小

我的问题是试图实现用于梯形总和其为下式Δx/2[y0 + 2y1 + 2y2 + … + 2yn-1 + yn]

这里是我的代码，我被困在地区是“求和”的一部分FOR循环。我想总结一下2y2 + 2y3....2yn-1因为我已经占了2y1 。我得到一个答案，但它并不像它应该的那样精确。例如，我得到6.071717974723753代替6.101605982576467 。

谢谢你的帮助！

function t=trapintegral(f,a,b,tol)
format compact; format long;
syms x;
oldtrap = ((b-a)/2)*(f(a)+f(b));
n = 2;
h = (b-a)/n;
newtrap = (h/2)*(f(a)+(2*f(a+h))+f(b));
while (abs(newtrap-oldtrap)>=tol)
    oldtrap = newtrap;
    for i=[3:n]
        dx = (b-a)/n;
        trapezoidsum = (dx/2)*(f(x) + (2*sum(f(a+(3:n-1))))+f(b));
        newtrap = trapezoidsum;
    end
end
t = newtrap;
end

为什么这个代码不工作的原因是因为有您的总和为梯形规则二细微的误差。什么我正是指的是这样一句话：

trapezoidsum = (dx/2)*(f(x) + (2*sum(f(a+(3:n-1))))+f(b));

召回梯形积分规则公式：

^{资料来源：维基百科}

对于第一个错误， f(x)应当f(a)如要包括的出发点，而不是应该保留象征性的。事实上，你应该简单地摆脱了syms x声明，因为它是不是在你的脚本非常有用。 a对应于x1通过咨询上述等式。

接下来的错误是第二届。实际上，你需要乘以你的指数值(3:n-1)通过dx 。此外，这实际上应该从去(1:n-1)和我稍后会解释。上述公式2去N ，但我们的目的，我们会从1到去N-1你有你的代码中设置了这样。

请记住，在梯形规则，你是细分有限区间为n件。所述^第一块i定义为：

x_i = a + dx*i;  ,

其中i去从1到N-1 请注意，这从1开始，而不是3就是因为第一块已经被考虑的原因， f(a) ，我们只数到N-1作为一块N由占f(b) 对于公式，这正好从2至N并通过修改代码这种方式，这正是我们到底在做什么。

因此，您的发言实际上需要是：

trapezoidsum = (dx/2)*(f(a) + (2*sum(f(a+dx*(1:n-1))))+f(b));

试试这个，让我知道，如果你得到正确的答案。 FWIW，MATLAB已经做实现梯形积分trapz作为@ADonda已经指出。然而，你需要正确建构你的x和y值设置此之前。换句话说，你需要设置你的dx手之前，然后计算出你的x使用点x_i我上面规定的公式，然后用这些来生成y值。然后，使用trapz计算面积。换一种说法：

dx = (b-a) / n;
x = a + dx*(0:n);
y = f(x);
trapezoidsum = trapz(x,y);

您可以使用上面的代码作为参考，看看你是正确地实现梯形规则。你的实现和使用上面的代码应该产生相同的结果。所有您需要做的是改变的价值n ，然后运行该代码为您的曲线下的不同细分领域的逼近。

编辑 - 2014年8月17日

我想通了，为什么你的代码不能正常工作。下面是原因：

该for环是不必要的。看看在for循环迭代。你有一个循环从去i = [3:n]你却不引用i的变量都在你的循环。因此，你不需要这个的。
你是不是正确计算连续间隔。你需要做的是，当你计算的梯形总和n个子区间，你再增加的这个值n ，然后再计算梯形规则。此值未在适当的增加while循环，这就是为什么你的区域是永远不会提高。
您需要保存里面的前一次的区域while循环，那么当你计算下一个区域，当你决定了区域之间的差异是否不小于公差的。我们也可以摆脱这些代码在试图开始并计算面积n = 2 。这是没有必要的，因为我们可以把这个里面你while循环。因此，这是你的代码应该是什么样子：

function t=trapintegral(f,a,b,tol)
format long; %// Got rid of format compact. Useless
%// n starts at 2 - Also removed syms x - Useless statement
n = 2;
newtrap = ((b-a)/2)*(f(a) + f(b)); %// Initialize
oldtrap = 0; %// Initialize to 0
while (abs(newtrap-oldtrap)>=tol)
    oldtrap = newtrap; %//Save the old area from the previous iteration
    dx = (b-a)/n; %//Compute width
    %//Determine sum
    trapezoidsum = (dx/2)*(f(a) + (2*sum(f(a+dx*(1:n-1))))+f(b));
    newtrap = trapezoidsum; % //This is the new sum
    n = n + 1; % //Go to the next value of n
end
t = newtrap;
end

通过运行代码，这是我得到：

trapezoidsum = trapintegral(@(x) (x+x.^2).^(1/3),1,4,0.00001)

trapezoidsum = 

6.111776299189033

警告

看看我定义你的函数的方式。您必须使用元素的元素作战为sum内环路命令将被量化。看看的^操作具体。您需要预先设置一个点来操作。一旦你这样做，我得到正确的答案。

编辑＃2 - 2014年8月18日

你说你要至少一个for循环。这是非常低效的，并且指定谁具有一种for代码回路真的不知道MATLAB是如何工作的。不过，你可以使用for循环积累的sum项。因此：

function t=trapintegral(f,a,b,tol)
format long; %// Got rid of format compact. Useless
%// n starts at 3 - Also removed syms x - Useless statement
n = 3;
%// Compute for n = 2 first, then proceed if we don't get a better
%// difference tolerance
newtrap = ((b-a)/2)*(f(a) + f(b)); %// Initialize
oldtrap = 0; %// Initialize to 0
while (abs(newtrap-oldtrap)>=tol)
    oldtrap = newtrap; %//Save the old area from the previous iteration
    dx = (b-a)/n; %//Compute width
    %//Determine sum
    %// Initialize
    trapezoidsum = (dx/2)*(f(a) + f(b));

    %// Accumulate sum terms
    %// Note that we multiply each term by (dx/2), but because of the
    %// factor of 2 for each of these terms, these cancel and we thus have dx
    for n2 = 1 : n-1
        trapezoidsum = trapezoidsum + dx*f(a + dx*n2);
    end
    newtrap = trapezoidsum; % //This is the new sum
    n = n + 1; % //Go to the next value of n
end
t = newtrap;
end

祝好运！