我的辛普森算法有什么问题?

Djamillah-虽然初始化了 h 可能仅适用于以下情况 a==0 ,因此您可能希望将这行代码更改为

h = (b-a)/(2^n+1);

我想知道 x = a:h:b; 将始终有效-有时 b 可能包含在列表中,有时可能不包含,具体取决于 小时 。您可能需要重新考虑并使用 linspace 相反

x = linspace(a,b,2^n+1);

这将保证 x 在间隔中均匀分布有2^n+1个点 [a,b] . 小时 然后可以初始化为

h = x(2)-x(1);

此外,在确定偶数和奇数索引时,我们需要忽略 x(x) 对于偶数和奇数。所以不是

xodd  = x(2:2:end-1);
xeven = x(3:2:end);

做

xodd  = x(2:2:end-1);
xeven = x(3:2:end-1);

最后,而不是使用向量 y (这是如何设置的?)我可能只是使用函数句柄来代替我正在集成的函数,并将上面的计算替换为

Isimp = delta/3*(func(x(1)) + 4*sum(func(xodd)) + 2*sum(func(xeven)) + ...
        func(x(end)));

除了这些微小的东西(可能可以忽略不计)之外,您的算法中没有任何东西表明存在问题。它产生了与我的版本类似的结果。

至于收敛的顺序 O(n^4) 或 O(h^4) ?

考虑到杰夫的建议,再做一些其他的改变,一切都如预期的那样。

%Reference solution with Simpson's method (the reference solution works well)
a=0;
b=1;
y=@(x) cos(x);
nmax=10;

%Simpson's method
Ehmatrix = [];
for n = 0:nmax
    x = linspace(a,b,2^n+1);
    h = x(2)-x(1);
    xodd = x(2:2:end-1);
    xeven = x(3:2:end-1);
    yodd = y(xodd);
    yeven = y(xeven);
    Isimp = (h/3)*(y(x(1))+4*sum(yodd)+2*sum(yeven)+y(b));
    E = abs(Isimp-integral(y,0,1));
    Ehmatrix(n+1,:) = [E h];
end

loglog(Ehmatrix(:,2),Ehmatrix(:,1))

P=polyfit(log(Ehmatrix(:,2)),log(Ehmatrix(:,1)),1);    
OrderofAccuracy=P(1)

你得到了 O(h) 准确性,因为 xeven=x(3:2:end) 是错误的。替换为 xeven=x(3:e:end-1) 修复代码,从而提高准确性。