複合辛普森規則的MatLab算法

我曾嘗試過，爲了它的樂趣，爲複合辛普森規則編寫一個MatLab代碼。據我所知，代碼是正確的，但我的答案並不像我想的那麼準確。如果我在函數f = cos（x）+ e ^（x^2）上嘗試我的代碼，其中a = 0，b = 1和n = 7，我的答案大約是1.9，當它應該是2時， 3。如果我使用維基百科可用的算法，我會得到一個非常接近的n = 7的近似值，所以我的代碼顯然不夠好。如果有人可以在我的代碼中看到任何錯誤，我會非常感激！複合辛普森規則的MatLab算法

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

來源

2012-10-29 Kristian

你的間隔[a,b]應分成n區間。這導致x的n+1值形成每個分區的邊界。您的向量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-29 21:24:02 mathematician1975

謝謝！我會研究這一點。 MatLab不會以x（0）開始索引的事實決不會永遠不會讓我惱火！ – Kristian

@Kristian是的，它可以讓生活變得困惑。看起來你至少需要根據公式使'x'包含'n + 1'元素，並且'x（0）= a'不是'f（a）'。同樣'x（n）= b'不'f（b）'（使用基於零的索引!!） – mathematician1975

再次感謝。我已經嘗試根據您的輸入修復我的代碼，但我的答案仍然有問題，恐怕。我非常確定索引在這裏拋出我，因爲我將我的代碼放在一個爲x（0）處的index-start編寫的僞代碼上。 – Kristian