C语言 用六种方法求定积分

C语言实验报告

求定积分

班级10信息与计算科学一班姓名戴良伟

学号2010750221

1. 描述问题

利用①左矩形公式，②中矩形公式，③右矩形公式 ，④梯形公式，⑤simpson 公式，⑥Gauss 积分公式求解定积分。

2. 分析问题

2.1定积分

21.1定积分的定义

定积分就是求函数()f x 在区间[],a b 中图线下包围的面积。即

()0,,,y x a x b y f x ====所包围的面积。这个图形称为曲边梯形，特例是曲边梯形。如下图：

（图1）

设一元函数()y f x =,在区间[],a b 内有定义。将区间[],a b 分成n 个小区间[][][][]00112,,,,,......,i a x x x x x x b 。设1i i i x x x -?=-，取区间i x ?中曲线上任意一点记做()i f ξ，作和式：

()1lim n n i f i xi ξ→+∞

=??

? ???

∑ 若记λ为这些小区间中的最长者。当0λ→时，若此和式的极限存在，则称这个和式是函数()f x 在区间[],a b 上的定积分。

记作：()b

a f x dx ?

其中称a 为积分下限，b 为积分上限，()f x 为被积函数，()f x dx 为被积式，∫ 为积分号。

之所以称其为定积分，是因为它积分后得出的值是确定的，是一个数，而不是一个函数。 21.2定积分的几何意义[1]

它是介于x 轴、函数f(x)的图形及两条直线x=a ，x=b 之间的各个部分面积的代数和。在x 轴上方的面积取正号；在x 轴下方的面积取负号。如图

2.2言实现定积分计算的算法

22.1利用复合梯形公式实现定积分的计算

假设被积函数为()f x ，积分区间为[],a b ，把区间[],a b 等分成n 个小区间，各个区间的长度为h ，即()/h b a n =-，称之为“步长”。根据定积分的定义及几何意义，定积分就是求函数()f x 在区间[],a b 中图线下包围的面积。将积分区间n 等分，各子区间的面积近似等于梯形的面积，面积的计算运用梯形公式求解，再累加各区间的面积，所得的和近似等于被积函数的积分值，n 越大，所得结果越精确。以上就是利用复合梯形公式实现定积分的计算的算法思想。

复合梯形公式：

()()()1

122n n i i h T f a f x f b -=??=++ ???

∑[2]

具体算法如下：

算法一

1：输入积分区间的端点值a 和b ；

2：输入区间的等分个数n （要求n 尽可能大，以保证程序运行结果有较高的精确度）；

3：计算步长()/h b a n =-；

4：对累加和赋初值()/2a b T f f =-； 5：计算累加和

()1

1n i i T f x -==∑

6：算出积分值n T T h =?； 7：输出积分近似值n T ，完毕。

1.2.2利用Smpson 公式实现定积分的计算

假设被积函数为()f x ，积分区间为[],a b ，把区间[],a b 等分成n 个小区间，各个区间的长度为h 。在复合梯形公式的基础上，构造出一种加速计算积分的方法。作为一种外推算法， 它在不增加计算量的前提下提高了误差的精度。 具体算法如下：

算法二

1：输入积分上限b 和下限a ；

2：输入区间的等分个数n （要求n 尽可能大，以保证程序运行结果有较高的精确度）；

3：利用辛甫生公式：[][][]()42/3S n T n T n =?-[2]，实现对定积分的求解（其中

[]2

T n均为梯形公式计算所得的结果，由此可见辛甫生公式是以梯形公式T n，[]

为基础的）；

；

4：算出积分值S

n

5：输出积分近似值S

，完毕。

n

1.2.3利用Guass公式实现定积分计算

Guass型求积公式是构造高精度差值积分的最好方法之一。他是通过让节点和积分系数待定让函数f(x)以此取i=0,1,2....n次多项式使其尽可能多的能够精确成立来求出积分节点和积分系数。高斯积分的代数精度是2n-1，而且是最高的。通常运用的是-1---+1的积分节点和积分系数，其他积分域是通过变换x=(b-a)t/2 +(a+b)/2 变换到-1到1之间积分。

算法三

1：输入积分上限b和下限a；

2：利用Guass公式，求定积分

4：算出积分值S

；

n

，完毕。

5：输出积分近似值S

n

3.程序的编写

3.1程序一（左矩形公式）

3.1.1源程序

#include

#include

void main()

{double f(double x);

/*f(x)为函数举例，即被积函数*/

int i,n;

/*n为区间等分的个数，应尽可能大*/

double a,b,h,s;

/*a为积分下限，b为积分上限，h为步长*/

printf("积分下限a:\n");

scanf("%lf",&a);

printf("积分上限b:\n");

scanf("%lf",&b);

printf("区间等分个数n :\n");

scanf("%d",&n);

h=(b-a)/n; /*步长的计算*/

s=f(a)*h;

for(i=1;i

{s=s+f(a+i*h)*h;

}

printf("函数f(x) 的积分值为s=%10.6f\n",s); }

/*以下为被积函数的定义，即函数举例*/ double f(double x)

{double y;

y=sqrt(4-x*x);

return (y);}

3.1.2程序一的编译运行

被积函数为f(x)=sqrt4-(x*x)的情况

先编译，再运行，屏幕显示及操作如下：

输入0+回车

输入2+回车

输入1000+回车

3.2程序二（中矩形公式）

3.2.1源程序

#include

#include

void main()

{double f(double x);

/*f(x)为函数举例，即被积函数*/

int i,n;

/*n为区间等分的个数，应尽可能大*/ double a,b,h,s;

/*a为积分下限，b为积分上限，h为步长*/ printf("积分下限a:\n");

scanf("%lf",&a);

printf("积分上限b:\n");

scanf("%lf",&b);

printf("区间等分个数n :\n");

scanf("%d",&n);

h=(b-a)/n; /*步长的计算*/

s=0.5*(f(a)+f(a+h))*h;

for(i=1;i

{s=s+0.5*(f(a+i*h)+f(a+(i+1)*h))*h;

}

printf("函数f(x) 的积分值为s=%10.6f\n",s); }

/*以下为被积函数的定义，即函数举例*/ double f(double x)

{double y;

y=sqrt(4-x*x);

return (y);}

3.2.2程序二的编译运行

被积函数为f(x)=sqrt4-(x*x)的情况

先编译，再运行，屏幕显示及操作如下：

输入0+回车

输入2+回车

输入1000+回车

3.3程序三（右矩形公式）

3.3.1源程序

#include

#include

void main()

{double f(double x);

/*f(x)为函数举例，即被积函数*/

int i,n;

/*n为区间等分的个数，应尽可能大*/ double a,b,h,s;

/*a为积分下限，b为积分上限，h为步长*/ printf("积分下限a:\n");

scanf("%lf",&a);

printf("积分上限b:\n");

scanf("%lf",&b);

printf("区间等分个数n :\n");

scanf("%d",&n);

h=(b-a)/n; /*步长的计算*/

s=f(a+h)*h;

for(i=1;i

{s=s+f(a+(i+1*h))*h;

}

printf("函数f(x) 的积分值为s=%10.6f\n",s); }

/*以下为被积函数的定义，即函数举例*/ double f(double x)

{double y;

y=sqrt(4-x*x);

return (y);}

3.3.2程序三的编译运行

被积函数为f(x)=sqrt4-(x*x)的情况

先编译，再运行，屏幕显示及操作如下：

输入0+回车

输入2+回车

输入1000+回车

3.4程序四（梯形公式）

3.4.1源程序

#include

#include

void main()

{double f(double x);

/*f(x)为函数举例，即被积函数*/

int i,n;

/*n为区间等分的个数，应尽可能大*/ double a,b,h,s;

/*a为积分下限，b为积分上限，h为步长*/ printf("积分下限a:\n");

scanf("%lf",&a);

printf("积分上限b:\n");

scanf("%lf",&b);

printf("区间等分个数n :\n");

scanf("%d",&n);

h=(b-a)/n; /*步长的计算*/

s=0.5*(f(a)+f(a+h))*h;

for(i=1;i

{s=s+0.5*(f(a+i*h)+f(a+(i+1)*h))*h;

}

printf("函数f(x) 的积分值为s=%10.6f\n",s); }

/*以下为被积函数的定义，即函数举例*/ double f(double x)

{double y;

y=sqrt(4-x*x);

return (y);}

3.4.2程序四的编译运行

被积函数为f(x)=sqrt4-(x*x)的情况

先编译，再运行，屏幕显示及操作如下：

输入0+回车

输入2+回车

输入1000+回车

3.5程序五（Simpson公式）

3.5.1源程序

#include

#include

void main()

{double T(double x,double y,int z); double a,b,s;

int n;

printf("积分下限a:\n");

scanf("%lf",&a);

printf("积分上限b:\n");

scanf("%lf",&b);

printf("区间等分个数n :\n");

scanf("%d",&n);

s=(4*T(a,b,2*n)-T(a,b,n))/3;

/*利用辛甫生公式求解定积分*/

printf("函数f(x) 的积分值为s=%f\n",s); }

/*以下为复合梯形公式的定义*/

double T(double x,double y,int z) {double h,Tn;

int i;

double f(double t);

h=(y-x)/z;

Tn=(f(x)+f(y))/2;

for(i=1;i

Tn=Tn+f(x+i*h);

Tn=Tn*h;

return (Tn);

}

/*以下为被积函数的定义，即函数举例*/ double f(double t)

{double s;

s=sqrt(4-t*t);

return (s);}

3.5.2程序四的编译运行

被积函数为f(x)=sqrt4-(x*x)的情况

先编译，再运行，屏幕显示及操作如下：

输入0+回车

输入2+回车

输入1000+回车

3.6程序六（Guass公式）

3.6.1源程序

#include

#include

#define N 3

float gass_integral(float (*)(float),float,float,int);

void main()

{

float function_name(float);

float a,b;

printf("请输入积分上限b\n");

scanf("%f",&b);

printf("请输入积分下限a\n");

scanf("%f",&a);

float ans;

ans=gass_integral(function_name,a,b,N);

printf("ans=%f",ans);

}

//高斯求积：代数精度为2n-1. -1---+1 之间

float gass_integral(float (*func)(float x), float a,float b ,int n ) {

//高斯点及其求积系数列表-------------------------------------------------------------------------------------------

float x1[1]={0.0}; float A1[1]={2};

float x2[2]={-0.5573503,0.5573503}; float A2[2]={1,1};

float x3[3]={-0.7745967,0.0,0.7745967}; float A3[3]={0.555556,0.888889,0.555556};

float x4[4]={0.3399810,-0.3399810,0.8611363,-0.8611363};

float A4[4]={0.6521452,0.6521452,0.3478548,0.3478548}; float x5[5]={0.0,0.5384693,-0.5384693,0.9061799,-0.9061799}; float

A5[5]={0.5688889,0.4786287,0.4786287,0.2369269,0.2369269};

//----------------------------------------------------------------------------------------------------------------------------------

float * p,* t;

switch ( n)

{case 1 : p=x1;t=A1;break;

case 2 : p=x2;t=A2;break;

case 3 : p=x3;t=A3;break;

case 4 : p=x4;t=A4;break;

case 5 : p=x5;t=A5;break;

default : printf("intput wrong!");

}

float g=0;

for(int i=0;i

{g+=(*func)((b-a)*p[i]/2+(a+b)/2)*t[i];}

g*=(b-a)/2;

return g;

}

float function_name(float x)

{return (sqrt(4-x*x));

}

3.6.2程序四的编译运行

被积函数为f(x)=sqrt4-(x*x)的情况

先编译，再运行，屏幕显示及操作如下：

输入2+回车

输入0+回车

4误差分析

手工计算结果为：3.156173.，左矩形公式误差：0.39%，中矩形公式误差：0.46%，右矩形公式误差：0.52%，梯形公式误差：0.46%，辛普森公式和高斯公式误差几乎等于0，六个程序运行结果对比，在计算相同的函数f(x)=sqrt(4-x*x)的定积分，Simpson公式和Guass公式比矩形和梯形公式更可行，更有效。