题目描述

对某函数进行近似积分求解。

自适应辛普森积分 递归求解

#include<bits/stdc++.h>
using namespace std ;

const double  eps = 1e-12 ;

double f(double x)
{
    return sin(x)/x ;
}

double simpson(double l , double r)
{
    double mid = (l+r) / 2 ;
    return (r-l) * (f(l) + 4*f(mid) + f(r)) / 6 ;

}

double asr(double l , double r , double area)
{
    double mid = (l + r ) / 2 ;

    double left = simpson(l , mid) , right = simpson(mid, r) ;

    if(fabs(left + right - area) < eps ) return left + right ;
    else return asr( l ,mid , left) + asr(mid , r , right ) ; 
} 

int main() 
{
    double l ,r ;
    cin >> l >> r ;
    printf("%lf\n" , asr( l , r , simpson(l,r))) ;
    return 0 ;
}

扩展

复合梯形公式 不推荐，因为余项是二次，收敛较慢

#include<bits/stdc++.h>
using namespace std ;

const double  eps = 1e-12 ;

double f(double x)
{
    return sin(x)/x ;
}

double Trapezoid(double l , double r) 
{
    return (f(l) + f(r)) * (r-l) / 2 ; 
}

int main() 
{
    double l ,r ;
    scanf("%lf%lf",&l,&r) ;
    double h = (r-l)/200000 ;  // 迭代次数很高
    double res = 0 ;
    double last = l ; 
    for(int i = 0 ; i < 200000 ; i ++  ) 
    {
        res += Trapezoid(last , last + h) ;
        last += h ;
    }
    printf("%lf\n" , res ) ;
    return 0 ;
}

复合辛普森积分 有着广泛应用，余项是四次，收敛快，编码简单，推荐

将广区间分割若干份，对每一份进行一次辛普森积分。

优势：去除一个递归函数，简化编码。
注意：步长h要足够小，才可以保证万无一失。

#include<bits/stdc++.h>
using namespace std ;

const double  eps = 1e-12 ;

double f(double x)
{
    return sin(x)/x ;
}

double simpson(double l , double r)
{
    double mid = (l+r) / 2 ;
    return (r-l) * (f(l) + 4*f(mid) + f(r)) / 6 ;

}

int main() 
{
    double l ,r ;
    cin >> l >> r ;
    double h = (r-l)/350 ;
    double res = 0 ;
    double last = l ; 
    for(int i = 0 ; i <350 ; i ++  ) 
    {
        res += simpson(last , last + h) ;
        last += h ;
    }
    printf("%lf\n" , res ) ;
    return 0 ;
}

另外，类似的算法还有n=4时的科茨公式，详见牛顿-科茨公式

复合科茨公式 余项是六次，收敛比以上都快，但是公式相比不如辛普森好记忆

复合的思路同复合辛普森，只是换了公式。

#include<bits/stdc++.h>
using namespace std ;

const double  eps = 1e-12 ;

double f(double x)
{
    return sin(x)/x ;
}

double Coates(double l , double r) 
{
    double d = (r-l)/4 ;
    double x0 = l ,x1 = l + d , x2 = x1+ d ,x3 = r -  d , x4 = r;
    return (r-l)*(7*f(x0) + 32*f(x1) +12*f(x2) +32*f(x3) +7*f(x4) ) / 90  ; 
}

int main() 
{
    double l ,r ;
    scanf("%lf%lf",&l,&r) ;
    double h = (r-l)/100 ;
    double res = 0 ;
    double last = l ; 
    for(int i = 0 ; i < 100 ; i ++  )     // 迭代次数极大降低 ，最低83次迭代
    {
        res += Coates(last , last + h) ;
        last += h ;
    }
    printf("%lf\n" , res ) ;
    return 0 ;
}

龙贝格公式 收敛快，最低87次迭代 （58 + 29）

利用科茨公式的外推形式得到

#include<bits/stdc++.h>
using namespace std ;

const double  eps = 1e-12 ;

double f(double x)
{
    return sin(x)/x ;
}

double Coates(double l , double r) 
{
    double d = (r-l)/4 ;
    double x0 = l ,x1 = l + d , x2 = x1+ d ,x3 = r -  d , x4 = r;
    return (r-l)*(7*f(x0) + 32*f(x1) +12*f(x2) +32*f(x3) +7*f(x4) ) / 90  ; 
}

int main() 
{
    double l ,r ;
    scanf("%lf%lf",&l,&r) ;
    double h = (r-l)/60 ;
    double res1 = 0 ;
    double last = l ; 
    for(int i = 0 ; i < 60 ; i ++  ) 
    {
        res1 += Coates(last , last + h) ;
        last += h ;
    }

    double res2 = 0 ;
    last = l ;
    h =  (r-l)/30 ;
    for(int i = 0 ; i < 30; i ++  ) 
    {
        res2 += Coates(last , last + h) ;
        last += h ;
    }

    printf("%lf\n" , res1 * 64 / 63  - res2 / 63  ) ;
    return 0 ;
}

总结：使用复合求法的话，如果记得住科茨公式，再搭配上比较高次数的迭代（比如1000次），基本万无一失了