此代码不起作用的原因是您对梯形规则的求和中有两个小错误。我所指的正是这种说法：
trapezoidsum = (dx/2)*(f(x) + (2*sum(f(a+(3:n-1))))+f(b));
回忆梯形积分规则的方程：
来源：Wikipedia
对于第一个错误，f(x) 应该是 f(a)，因为您包含了起点，并且不应保留为符号。事实上，您应该简单地去掉 syms x 语句，因为它在您的脚本中没有用处。 a 对应于x1，通过查阅上述等式。
下一个错误是第二项。您实际上需要将索引值(3:n-1) 乘以dx。此外，这实际上应该来自(1:n-1)，我稍后会解释。上面的等式从 2 变为 N，但出于我们的目的，我们将从 1 变为 N-1，因为您的代码设置是这样的。
请记住，在梯形规则中，您将有限区间细分为n 部分。第 ith 块定义为：
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 已经通过 @ADonda 已经指出的 trapz 实现了梯形集成。但是，在设置之前，您需要正确构建 x 和 y 值。换句话说，您需要事先设置您的dx，然后使用我上面指定的x_i 等式计算您的x 点，然后使用这些来生成您的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 变量。因此，您根本不需要这个。
您没有正确计算连续间隔。您需要做的是在计算nth 子区间的梯形和时，然后增加 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
祝你好运！