※ 引述《LESPAUL (大帥)》之銘言：
: 請教各位高手
: 如何不使用內建的積分quad積分指令來寫一個辛浦森積分
: 我想寫 W在[a,b]區間之內，z及Wn都為一個固定常數下
: 對(z.*(W/Wn).^2)/((1-(W/Wn).^2).^2+(2*z.*(W/Wn)).^2)做積分
: 如何使用辛浦森積分公式
: b
: S s(x) dx= h/3*{s + 4(s +s +s +...+s ) + 2(s +s +...s )+s }
: a 0 1 3 5 2n-1 2 4 2n-2 2n
: 其中a,b分別表示積分的下限與上限
: n指ab間成對的區間數目
: h=(b-a)/(2n) 代表每區間的半寬度
: 請教各位高手要怎麼寫辛浦森的積分
: 感激不盡啊
新浦森法有一個限制
就是他拿來積分的點個數一定要是奇數個（觀察一下你方程式的index就可以知道）
也就是使用時取的區間必須為偶數個
所以只要注意一下index的特色就很好寫了
s = 5 % 任意值，你想切幾等分一半的數字就是
N = 2 * s ; % 只是確定區間為偶數
xn = linspace( a , b , N + 1 ) ;
yn = f( xn ) ; % 你要的方程式給xn輸出成yn，記得yn要輸出成矩陣
Is = (xn(2)-xn(1))*(yn(1)+yn(end)+4*sum(yn(2:2:end))+2*sum(yn(3:2:end-1)))/3 ;
這樣積分就出來了
這只是最簡單的給幾個區間就直接計算出答案的寫法
如果要進行收斂測試的話
區間數量的變化必須為2 4 8 16 32 64 ......
這邊也只需要小調整即可
不過要做到收斂測試的話辛浦森法其實比較沒啥作用了
等於只是Romberg method的第二層而已
而第二層卻可以使用梯形法計算的結果處理出來
也就是說
想要玩收斂測試的話只需要梯形法配合Romberg method
根本不必再寫一個新浦森法
--
Deserves death! I daresay he does. Many that live deserve death. And some die
that deserve life. Can you give that to them? Then be not too eager to deal out
death in the name of justice, fearing for your own safty. Even the wise cannot
see all ends.
Gandalf to Frodo
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 140.120.25.235
※ 編輯: Gwaewluin 來自: 140.120.25.235 (11/21 10:53)