1. 程式人生 > >Romberg法求定積分

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

請按任意鍵繼續. . .