复合辛普森规则的 MatLab 算法答案

【问题标题】：MatLab algorithm for composite Simpson's rule复合辛普森规则的 MatLab 算法
【发布时间】：2012-10-19 06:00:13
【问题描述】：

为了好玩，我已经尝试为复合辛普森规则编写 MatLab 代码。据我所知，代码是正确的，但我的答案并不像我想要的那样准确。如果我在函数 f = cos(x) + e^(x^2) 上尝试我的代码，a = 0，b = 1 和 n = 7，我的答案大约是 1,9，而它应该是 2， 3.如果我使用 Wikipedia 上提供的算法，我会得到一个非常接近 n = 7 的近似值，所以我的代码显然不够好。如果有人能看到我的代码中的任何错误，我将不胜感激！

function x = compsimp(a,b,n,f)
% The function implements the composite Simpson's rule

h = (b-a)/n;
x = zeros(1,n+1);
x(1) = a;
x(n+1) = b;
p = 0;
q = 0;

% Define the x-vector
for i = 2:n
    x(i) = a + (i-1)*h;
end

% Define the terms to be multiplied by 4
for i = 2:((n+1)/2)
    p = p + (f(x(2*i -2)));
end

% Define the terms to be multiplied by 2
for i = 2:((n-1)/2)
    q = q + (f(x(2*i -1)));
end

% Calculate final output
x = (h/3)*(f(a) + 2*q + 4*p + f(b));

【问题讨论】：

标签： matlab numerical-integration

【解决方案1】：

更正应在p1、p2 和p3 部分。我试了一下，得到了大致准确的结果：

【讨论】：

我喜欢“近似精确”这个词... :P

【解决方案2】：

function x = compsimp(a,b,n,f)

我不知道这是否重要但不应该是f的第一个字母：

function x = compsimp(f,a,b,n)

【讨论】：

【解决方案3】：

您创建的代码运行良好。我看到的唯一问题是 n。根据我的经验，任何函数都可以尝试 n>=10000。

【讨论】：

复合辛普森规则所犯的错误（绝对值）由 \tfrac{h^4}{180}(ba) \max_{\xi\in[a,b]} |f^{(4)}(\xi)|，其中 h 是“步长”，由 h=(ba)/n 给出。
你应该edit你的答案而不是发布cmets。评论更多是为了与其他 SO 成员进行澄清/讨论。

【解决方案4】：

您的区间 [a,b] 应拆分为 n 区间。这导致形成每个分区边界的x 的n+1 值。您的向量 x 仅包含 n 元素。您的代码似乎只根据需要处理 n 术语而不是 n+1。

编辑:: 现在你已经根据上面的修改了问题，试试这个

% The 2 factor terms
for i = 2:(((n+1)/2) - 1 )
    q = q + (f(x(2*i)));
end

% The 4 factor terms
for i = 2:((n+1)/2)
   p = p + (f(x(2*i -1)));
end

【讨论】：

谢谢！我会调查这个。 MatLab 没有从 x(0) 开始索引的事实永远不会让我感到厌烦！
@Kristian 是的，它会让生活变得混乱。看来至少您需要根据公式使x 包含n+1 元素，还有x(0) = a 而不是f(a)。同样x(n) = b 而不是f(b)（使用基于零的索引！！）
再次感谢。我已经尝试根据您的输入修复我的代码，但恐怕我的答案还有很长的路要走。我很确定索引将我甩在这里，因为我的代码基于为 x(0) 处的 index-start 编写的伪代码。
@Kristian 为什么不将 x 划分为 n 个区间，然后使用 x 上的循环交替添加在内部 x 点处评估的函数的 2 和 4 的倍数，而不是分别计算每个总和。这可能会减少一些索引混乱。或者使用基于 1 的索引而不是基于零的索引写出辛普森一家公式，以便您的代码和公式索引匹配？？
@Kristian 我更新的答案有帮助吗？我尝试根据 wiki 公式基于 1 的索引重新计算公式