0

楽しみのために、複合シンプソンの法則のMatLabコードを書いてみました。私が見る限り、コードは正しいですが、私の答えは私が望むほど正確ではありません。関数f=cos（x）+ e ^（x ^ 2）でコードを試してみると、a = 0、b = 1、n = 7で、答えはおよそ1,9ですが、2になるはずです。 3.3。ウィキペディアで利用可能なアルゴリズムを使用すると、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 
function x = compsimp(a,b,n,f)

これが重要かどうかはわかりませんがf、最初の文字にするべきではありません。

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

修正は、、、およびの部分p1で行う必要があります。私はそれを試し、ほぼ正確な結果を得ました：p2p3

ここに画像の説明を入力してください

于 2016-11-03T01:25:59.190 に答える