代码之家 › 专栏 › 技术社区 › Djamillah

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

numerical-methods matlab

Djamillah · 技术社区 · 11 年前

我试图用Simpson的方法写一个近似积分的算法。然而,当我试图将其绘制成对数图时,我没有得到正确的精度顺序,即O(h^4)(我得到O(n))。但我找不到任何错误。这是我的代码:

%Reference solution with Simpson's method (the reference solution works well)
yk = 0;
yj = 0;
href = 0.0001;
mref = (b-a)/href;

for k=2:2:mref-1
    yk=yk+y(k*href+a);
end

for j=1:2:mref
    yj=yj+y(href*j+a);
end
Iref=(href/3)*(y(a)+y(b)+2*yk+4*yj);

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

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

a和b是积分极限,y是我们要近似的被积函数。

2 回复 | 直到 11 年前

Geoff 11 年前

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) ?

David 11 年前

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

%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) 修复代码,从而提高准确性。