龍貝格求解數值積分
阿新 • • 發佈:2019-01-31
龍貝格求積公式,是一種梯形法的遞推化,也是將求積區間逐次二分,將前一次二分的積分和Tn與下一次二分的積分和T2n有個遞推關係,通過梯形公式的餘項泰勒展開成級數形式,變數代換,進行整體的加減組合,能夠得到關於T2n和Tn的線性組合S。同理,對S的線性組合可以推出C,對C的線性組合,可以推出R。精度在逐次提高,每一次遞推提高2階。
下面是實現程式碼:
/*龍貝格求積分
*double romb(double a,double b,double eps,double(*f)(double))
*double a:給出積分割槽間下限
*double b:積分割槽間上限
*double eps:誤差精度
*double (*f)(double) :被積函式
*函式返回double
*/
#include<iostream>
#include<cmath>
double romb(double a, double b, double eps, double(*f)(double))
{
int i, k;
double y[10], p, x, q;
double h = b - a;
y[0] = h*((*f)(a) + (*f)(b)) / 2.0;
int m = 1,n=1;
double ep = eps + 1.0;
while ((ep >= eps) && (m <= 9 ))
{
p = 0.0;
for (i = 0; i < n; i++)
{
x = a + (i + 0.5)*h;
p = p + (*f)(x);
}
p = (y[0] + h*p) / 2.0;
for (k = 1; k <= m; k++)
{
q = (4.0*p - y[k - 1]) / 3.0;
y[k - 1] = p;
p = q;
}
ep = fabs (q - y[m - 1]);
m = m + 1;
y[m - 1] = q;
n = n * 2;
h = h / 2.0;
}
return q;
}
//積分函式:f(x)=x/(4+x*x)
double func(double x)
{
return x / (4 + x*x);
}
int main()
{
double a = 0.0;
double b = 1.0;
double eps = 0.0000001;
double t = romb(a, b, eps, func);
std::cout << "the value is " << t << std::endl;
return 0;
}