Tauruser

先介紹一下Romberg求積。
6.3　外推原理與Romberg求積

6.3.1　復合梯形公式遞推化與節點加密

　　在計算機上用等距節點求積公式時，若精度不夠可以逐步加密節點.設將區間分為n等分，節點，在區間上梯形公式為
　　　　　　　　　　　　　
若節點加密一倍，區間長為，記中點為在同一區間上的復合梯形公式惟
　　　　　　　　　　　　　
于是
　　　　　　　(6.3.1)
它表明是在的基礎上再加新節點的函數值之和乘新區間長，而不必用(6.2.6)重新計算，這時有誤差估計式
　　　　　　　　　　　　　
若，則得
　　　　　　　　　　　　(6.3.2)
它表明用，其誤差近似.這也是在計算機上估計梯形公式誤差的近似表達式.若(給定精度)，則.
　　若在區間［a,b］中做2n等分時，在上用Simpson公式計算，則由(6.2.8)可知
　　　
它恰好是(6.3.2)中I(f)的近似值，即

它表明用(6.3.2)計算I(f)，其精度已由提高到如果再將區間分半，使分為4個小區間，長度為，則可由(6.3.1)計算出及，利用復合公式余項(6.2.9)得
　　　　　　　　　
　　　　　　　　　
如果，則有
　　　　　　　　　　　　　　　(6.3.3)
從而有復合Simpson公式的誤差估計
　　　　　　　　　
如果用(6.3.3)近似，即
　　　　　　　　　　　　　　　(6.3.4)
則精度可達到.類似做法還可繼續下去.這樣對區間逐次分半，利用公式(6.3.1)逐次遞推.再由(6.3.2)，(6.3.3)逐次構造出精度愈來愈高的計算積分I(f)的公式，這就是Romberg求積的基本思想.

以下為我自己寫的求積程序。

// ?RombergIntegral.cpp?:?定義控制臺應用程序的入口點。
//

#include? < cmath >
#include? < iostream >
#include? < vector >
using ? namespace ?std;
const ? double ?PRECISION(. 000001 ); // 精度控制
const ?unsigned? int ?MAXK( 20 ); // 求解步驟控制
double ?RombergIntegral( double ?( * f)( double ?x), double ?a,? double ?b);
vector < vector < double >> ?T; // 用于存儲T表
double ?f( double ?x) // 要求的積分函數
{
???? return ?x * sin(x);
}
int ?_tmain( int ?argc,?_TCHAR * ?argv[])
{
????cout << " 本程序用于求解函數f(x)=x*sin(x)在0到6.28的積分 " << endl;
????cout << " 積分結果為： " << RombergIntegral(f, 0 , 6.28 ) << endl;
????cout << " 精度為 " << PRECISION << endl;
???? return ? 0 ;
}

double ?RombergIntegral( double ?( * f)( double ?x), double ?a,? double ?b)
{
???? int ?k( 0 );
???? double ?h = b - a;
????vector < double > ?temp;
????T.push_back(temp);
????T[ 0 ].push_back(h * (( * f)(a) + ( * f)(b)) / 2 );
???? for (k = 1 ; 1 ; ++ k)
???? {
????????T.push_back(temp);
????????T[ 0 ].push_back( 0.5 * T[ 0 ][k - 1 ]);
???????? for ( int ?i = 0 ;i < pow( 2 .,k - 1 ); ++ i)
???????? {
????????????T[ 0 ][k] += 0.5 * h * (( * f)(a + h / 2 + i * h));
????????}
???????? for ( int ?i = 1 ;i <= k; ++ i)
????????????T[i].push_back((pow( 4 .,i) * T[i - 1 ].back() - T[i - 1 ][T[i - 1 ].size() - 2 ]) / (pow( 4 .,i) - 1 ));
????????h /= 2 ;
???????? double ?temp = T[k].back();
???????? if (fabs(T[k].front() - T[k - 1 ].front()) < PRECISION? || ??k == MAXK)? break ; // 到
????}
????
???? return ?T[k].back();
}

// 以上程序在vs2005+win2003下編譯運行通過。