复合辛普森规则的MatLab算法

Question

我尝试过，只是为了好玩，为复合辛普森法则写了一个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));

Answer 1

您的间隔间隔应拆分为多个n [a,b]。这将导致形成每个分区边界的x的n+1值。向量x只包含n元素。看起来您的代码只处理n术语，而不是所需的n+1。

EDIT:：现在你已经根据上面的内容修改了问题，试试这个

% 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

Answer 2

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

Answer 3

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

我不知道这是否重要，但不应该在第一个字母中使用f：

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