数值积分-龙贝格(Romberg)积分

　　数值积分在工程上是个比较有用的数学工具。在工程上有很多数学问题，看似简单，计算所用的数学公式不算复杂，但是求解起来却很困难，很难获得解析解的公式，这个时候就需要用到数值求解的办法来获取满足工程需要的近似数值结果。

　　举个简单的例子，求一个变截面圆柱（或者叫圆台）的等效截面惯性矩，计算公式为Ie=H^3/(3∫[0,H](x^2/I(x))dx。这个公式使用的是结构力学里面的力法原理推导的，积分本来也不是很复杂，问题出在I(x)上，I(x)是沿柱子高度方向坐标x处截面的惯性矩,一般情况，I(x)是关于x的4次多项式。本来以为可以找到这个积分的解析表达式，试着积了一下，才发现很困难。去找积分公式，只找到2次的，1/(ax^2+bx+c)的积分，从结果来看，是要解方程的。4次方程的解析解是可以求出来，不过太复杂了。我试着用MatLab积分了一下，给出结果好长啊，有些符号我还不认识，不过代入数值还是有比较简单的结果。

 

syms x a b c d e Ie H f=x*x/(a*x^4+b*x^3+c*x^2+d*x+e) Ie=H^3/3/int(f,x,0,H) Ie = 1/3*H^3/(sum(_R^2/(4*a*_R^3+3*b*_R^2+2*c*_R+d)*log(H-_R), _R = RootOf(a*_Z^4+b*_Z^3+c*_Z^2+d*_Z+e))-sum(_R^2/(4*a*_R^3+3*b*_R^2+2*c*_R+d)*log(-_R), _R = RootOf(a*_Z^4+b*_Z^3+c*_Z^2+d*_Z+e)))  

　　最后无奈只好求助于数值积分来搞定这个小怪物。

//.h //龙贝格积分 double Romberg(double (*f)(double, CDoubleArray &), double x1, double x2, CDoubleArray & a); //.cpp double Romberg(double (*f)(double, CDoubleArray &), double x1, double x2, CDoubleArray & a) { //x1, x2 积分上下限 //a为f函数的参数表 int m, n;//m控制迭代次数, 而n控制复化梯形积分的分点数. n=2^m double h, x; double s, q; int MAXREPT = 10; double ZERO6 = 0.000001; double ep; //精度要求 CDoubleArray y; y.SetSize(MAXREPT); //每次循环依次存储Romberg计算表的每行元素,以供计算下一行,算完后更新 double p;//p总是指示待计算元素的前一个元素(同一行) //迭代初值 h = x2 - x1; y[0] = h*(f(x1, a) + f(x2, a))/2.0; m = 1; n = 1; ep = ZERO6 + 1.0; //迭代计算 while ((ep >= ZERO6) && (m < MAXREPT)) { //复化积分公式求T2n(Romberg计算表中的第一列),n初始为1,以后倍增 p = 0.0; for (int i=0; i<n; i++)//求Hn { x = x1 + (i+0.5)*h; p = p + f(x, a); } p = (y[0] + h*p)/2.0;//求T2n = 1/2(Tn+Hn),用p指示 //求第m行元素,根据Romberg计算表本行的前一个元素(p指示), //和上一行左上角元素(y[k-1]指示)求得. s = 1.0; for (int k=1; k<=m; k++) { s = 4.0*s; q = (s*p - y[k-1])/(s - 1.0); y[k-1] = p; p = q; } p = fabs(q - y[m-1]); m = m + 1; y[m-1] = q; n = n + n; h = h/2.0; } y.RemoveAll(); return q; }