0

我已经尝试过,只是为了好玩,为复合辛普森规则编写一个 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));
4

4 回答 4

1

你的区间[a,b]应该分成n区间。这导致形成每个分区边界的n+1值。x您的向量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
于 2012-10-29T21:24:02.503 回答
0

您创建的代码工作正常。我看到的唯一问题是 n。根据我的经验,尝试 n>=10000 用于任何功能。

于 2015-04-28T18:55:08.507 回答
0

更正的部分应该p1p2p3。我试了一下,得到了大致准确的结果:

在此处输入图像描述

于 2016-11-03T01:25:59.190 回答
0
function x = compsimp(a,b,n,f)

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

function x = compsimp(f,a,b,n)
于 2015-11-14T21:02:08.803 回答