Tauruser

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

6.3.1　復(fù)合梯形公式遞推化與節(jié)點(diǎn)加密

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

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

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

// ?RombergIntegral.cpp?:?定義控制臺應(yīng)用程序的入口點(diǎn)。
//

#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) // 要求的積分函數(shù)
{
???? return ?x * sin(x);
}
int ?_tmain( int ?argc,?_TCHAR * ?argv[])
{
????cout << " 本程序用于求解函數(shù)f(x)=x*sin(x)在0到6.28的積分 " << endl;
????cout << " 積分結(jié)果為： " << 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下編譯運(yùn)行通過。