对于积分:
运行结果如下图:
只要找到被积公式的原函数F(x),利用牛顿莱普利兹公式有:
但是,实际使用这种求积分的方法往往是有困难的,因为大量的被积函数的原函数是不能用初等函数表示的;另外,当f(x)是由测量或数值计算给出的一张数据表时,牛顿莱普利兹公式也无法直接运用,因此有必要研究积分的数值计算问题。
对于一些理论的推导,大家可以看看维基百科,下面我
主要给出牛顿-科特斯公式在n=1(梯形求积公式)、n=2(辛普森公式)的情况,并通过代码实现。
梯形公式:
辛普森公式:
应用高阶牛顿-科特斯公式计算积分时,会出现数值不稳定的情况,而低阶公式往往因为积分步长过大使得离散误差变大,因此,为了提高求积公式的精度,可以把积分区间分成若干个子区间,在每个子区间上使用低阶求积公式,然后将结果加起来,这种方法称为复化求积法。
复化梯形公式
将区间[a,b]划分为n等分,步长为h=(b-a)/h,节点为
,在每个子区间
使用梯形公式得:
复化辛普森公式
根据复化梯形公式的推导,同理可得复化辛普森公式为:
下面我们通过实例来实现复化梯形公式和复化辛普森公式:
对于函数f(x)=sin(x)/x,试用复化梯形公式和复化辛普森公式计算函数f(x)在[0,1]上的积分。
具体的程序实现如下:
#include<stdio.h> #include<math.h> double Function(double x)//所要计算积分的函数f(x) { if(x==0)//sin(x)/x在0处的取值为1 return 1; else return sin(x)/x; } //复化梯形公式 double Trapz(double a,double b,int n) { double h=(b-a)/n; double T=0; for(int i=1;i<n;i++) { T=T+Function(a+i*h); } T*=2; T=(Function(a)+Function(b)+T)*h/2; return T; } //复化辛普森公式 double MulripleSimpson(double a,double b,int n) { double h=(b-a)/n; double T=0; for(int i=0;i<n;i++) { T=T+Function(a+i*h)+4*Function(a+(i+0.5)*h)+Function(a+(i+1)*h); } T=T*h/6; return T; } void main() { printf("使用复化梯形公式可得:%f\n",Trapz(0,1,8)); printf("使用复化辛普森公式可得:%f\n",MulripleSimpson(0,1,4)); }
运行结果如下图:
结果分析:
比较复化梯形公式和复化辛普森公式两种方法的运行结果,我们发现复化辛普森公式与准确值0.9460831更加接近,复化梯形公式只有2位有效数字,而复化辛普森公式有6为有效数字。
原文:http://blog.csdn.net/tengweitw/article/details/43311685
作者:nineheadedbird