Romberg法求定積分
1 實驗題目：
2 演算法組織：
2.1 演算法原理：
由梯形演算法可以推得：
，這裡h為2n個區間的步長，即h = (b - a) / 2n;
如果有2^n個相同的子區間，則上式變為：
，這裡hn為2^n個區間的步長，即hn = (b - a) / 2^n;
所以，對於Romberg演算法，可以採用如下公式計算：
，
，
虛擬碼：
input a, b, M
h ← b - a
for n = 1 to M do
h ← h / 2
for m = 1 to n do
end do
end do
output R(n,m) (0 ≤ n ≤ M, 0 ≤ m ≤ n)
（演算法詳細說明可參考《數值分析（原書第3版）》，（美）David Kincaid, Ward Cheney著，王國榮、俞耀明、徐兆亮譯，機械工業出版社，2005年9月第1版，第400到402頁。）
2.2 程式碼:
程式設計環境為VS 2005。原始碼檔案為“上交用_7.cpp”。
依據2.1中的虛擬碼程式設計，並添加了誤差判斷條件，當|R[n][m]-R[n-1][m]| <= e(e為指定誤差)時，程式結束。
//ROMBERG求定積分
#include "stdafx.h"
#include <iostream>
#include <complex>
//exp(const complex& _ComplexNum)函式
#include //使用三角函式及log函式
const int M = 10;//用於指定ROMBERG函式中二維陣列的長度
double ROMBERG(double a, double b, double e);
double f(double x);
int main()
{
using namespace std;
cout << "Romberg積分結果為: " << ROMBERG(0, 1, 0.00001) << endl;
return 0;
}
double f(double x)
{
return (4.0 / (1 + x * x));
}
double ROMBERG(double a, double b, double e)
{
double R[M][M] = {0};
double h = b - a;//步長
double temp;//
R[0][0] = h * (f(a) + f(b)) / 2;
for(int n = 1; n < M; n++)
{
h /= 2;
temp = 0;
for(int i = 1; i <= pow(2.0,n-1); i++)
temp += f( a + ( 2 * i - 1 ) * h );
R[n][0] = R[n-1][0] / 2 + h * temp;
for(int m = 1; m < 4; m++)//只求到第四列
{
R[n][m] = R[n][m-1] + ( R[n][m-1] - R[n-1][m-1] ) / ( pow(4.0,m) - 1 );
//判斷誤差
if( n>=4 && m==3 )
if(fabs(R[n][m]-R[n-1][m]) <= e)
{
//輸出結果
for(int i = 0; i <= 3; i++)
{
for(int j = i; j <= n; j++)
std::cout << "R[" << j << "][" << i << "] = " << R[j][i] << std::endl;
std::cout << std ::endl;
}
return R[n][3];
}
}
}
return R[M][3];
}
3 計算結果
R[0][0] = 0.920735
R[1][0] = 0.939793
R[2][0] = 0.944514
R[3][0] = 0.945691
R[4][0] = 0.945985
R[1][1] = 0.946146
R[2][1] = 0.946087
R[3][1] = 0.946083
R[4][1] = 0.946083
R[2][2] = 0.946083
R[3][2] = 0.946083
R[4][2] = 0.946083
R[3][3] = 0.946083
R[4][3] = 0.946083
Romberg積分結果為: 0.946083
請按任意鍵繼續. . .