用梯形法求定积分的值

最新推荐文章于 2023-06-01 00:18:22 发布

rainbow3b

最新推荐文章于 2023-06-01 00:18:22 发布

阅读量6.7k

点赞数 8

分类专栏： C++算法学习文章标签：算法 c++

本文链接：https://blog.csdn.net/rainbow3b/article/details/128200479

版权

C++算法学习专栏收录该内容

6 篇文章 0 订阅

订阅专栏

一、梯形法求解定积分的过程

1.求定积分值存在的问题

计算定积分是数值计算领域内的一个重要内容。对于能够得到原函数的被积函数，如：

$\int xdx=\frac{1}{2}x^{2}$ ,

其定积分可以直接计算。

$\int_{0}^{1}xdx=\frac{1}{2}1^{2}-\frac{1}{2}0^{2}=\frac{1}{2}$

但对于不易得到原函数的被积函数，可以考虑使用数值计算的方法得到近似值。如：

$f(x)=\frac{sinx}{x}$

不易得到原函数，故其如下

$\int_{0}^{1}sinx/xdx$

的定积分也不容易求解。

2.定分积的几何意义

定积分的几何意义是被积函数和x轴以及积分上限、积分下限之间围成的图形的面积。如下图所示：

图中x轴、y=x与y=1所围三角形的面积即为

$\int_{0}^{1}xdx=\frac{1}{2}1^{2}-\frac{1}{2}0^{2}=\frac{1}{2}$

对于无法得到原函数的被积函数，其定积分也是这样的面积。如下图所示：

上图中，x轴、f(x)、a、b等所围成的阴影面积即为

$\int_{a}^{b}f(x)dx$

的值，即以a为积分下限、b为积分上限、被积函数f(x)的积分值。

3.积分数值计算过程

图2中的阴影面积是不容易得到的。可以用近似的方法得到，即梯形近似法。如下图所示：

在f(a)、f(b)之间直接连线，则可以由x轴、a、b、以及连线形成梯形。该梯形的两个底分别为：f(a)、f(b)，高为（b-a)，由梯形的面积公式，可以得到上图中所围梯形的面积为：

$S=\frac{(f(a)+f(b))*(b-a)}{2}$ (式1)

这个可以作为 $\int_{a}^{b}f(x)dx$ 的近似值。当然其误差比较大。

如果对上图进一步分割，按下图进行拆分，

设c为积分下限a、积分上限b的之间中点，则

$c=\frac{a+b}{2}$

以f(c)为线对图形进行分割，形成两个梯形，分别计算面积，再对面积求和。

梯形S1的两个底分别为f(a)、f(c)，高为(c-a)；梯形S2的两个底分别为f(c)、f(b)，高为(b-c)。则面积：

$S1=\frac{f(a)+f(c)}{2}*(c-a) =\frac{f(a)+f(c)}{2}*(\frac{b+a}{2}-a) =\frac{f(a)+f(c)}{2}*(\frac{b-a}{2}) =(\frac{b-a}{4})\({f(a)+f(c)}) =(\frac{b-a}{4})\({f(a)+f(\frac{b+a}{2})})$

同样，

$S2=(\frac{b-a}{4})\({f(b)+f(\frac{b+a}{2})})$

这样，总的面积是：

$S=S1+S2=(\frac{b-a}{4})\({f(a)+f(\frac{b+a}{2})})+(\frac{b-a}{4})\({f(b)+f(\frac{b+a}{2})})=(\frac{b-a}{4})\({f(a)+f(b)+2*f(\frac{b+a}{2})})$ (式2)

令（式1）为 $S_{k0}$ ，（式2）为 $S_{k1}$ ，可以得到：

$S_{k1}=\frac{1}{2}S_{k0}+\frac{b-a}{2}*f(\frac{b+a}{2})$ （式3）

如果对上图中S1、S2再进行等分分割，可以得到如下图的S1、S2、S3、S4四个梯形。再次分别计算其面积后进行求和，可以得到更近似的面积，即定积分值。如此，不断进行将[a,b]进行二分，递推下去，得到的积分值越来越精确，可以无限逼近函数f(x)在[a，b]上的定积分。

这样可以得到递推公式。具体是：将[a，b]段分为n等分，一共有n+1个等分点，第k(k=0,1,2...n)个等分点的值为： $x_{k}=a+k*h$ ，其中， $h=\frac{b-a}{n}$

用 $S_{n}$ 表示将[a,b]n等分后梯形面积的总各，用 $S_{2n}$ 表示将[a,b]2n等分后的梯形面积总各，得到递推公式为：

$S_{2n}=\frac{1}{2}*S_{n}+\frac{h}{2}*\sum_{k=0}^{n}f(\frac{1}{2}*(x_{k}+x_{k+1}))$ （式4）

利用上式，即可以得到定积分的近似值。设置某一个精度e，当 $\left | S_{2n}-S_{n} \right |<e$ 时，即认为达到精度，迭代结束。

二、计算定积分 $\int_{0}^{1}\frac{sinx}{x}dx$ 的C++程序

1.计算函数 $\frac{sinx}{x}$ 的值

//---------------------------
double func( double x )
{
	if( x != 0 )
		return sin( x ) / x;
	else
		return 1.0;
}
//----------------------------

由于x不能为0，且 $\lim_{x\rightarrow 0}\frac{sinx}{x}=1$ ，因此，设定，当x=0时，返回值为1.0。

2.迭代函数

//----------------------------
double ING( double a, double b, double e )//参数：积分下限，积分上限，精度
{
	double T1 = 0.0;
	double T2 = 0.0;
	double S = 0.0;
	double h, x;
	int flag;
	h = b - a;
	T1 = h / 2 * ( func( a ) + func( b ) );//计算第一个面积
	do {
		double s = 0;
		x = a + h /2;
		while( x < b )//循环计算
		{
			s = s + func( x );
			x = x + h;
		}
		T2 = T1 / 2 + h / 2 * s;
		if( fabs( T1 - T2 ) >= e )//判断精度
		{
			h = h / 2;
			T1 = T2;
			flag = 1;//当精度不够时，标志为1
		}
		else flag = 0;
	}
	while( flag );//为1时（意味着精度不够），继续计算。
	return T2;
}
//-------------------------------

3.完整程序

//-----------------------------
#include <iostream>
#include <cmath>
using namespace std;
//-----------------------------
double func( double x );
double ING( double a, double b, double e );
//-----------------------------
int main()
{
	double LowLim, HighLim, Accur;
	cout << "Input the low and high linitation and accuracy:" << endl;
	cout << "Low limitation: ";
	cin >> LowLim;           //输入积分下限
	cout << "High limitation:";
	cin >> HighLim;          //输入积分上限
	cout << " Accracy:";
	cin >> Accur;            //输入精度
	
	if( LowLim >= HighLim )  //积分上限不能小于积分下限，否则提醒再次输入
		cout << "HighLim is less than LowLim.Please input again.";
	else cout << "The result of integration is: " << ING( LowLim, HighLim, Accur ) << endl;
	system( "PAUSE" );
	return 0;
}
//---------------------------
double func( double x )
{
	if( x != 0 )
		return sin( x ) / x;
	else
		return 1.0;
}
//----------------------------
double ING( double a, double b, double e )//参数：积分下限，积分上限，精度
{
	double T1 = 0.0;
	double T2 = 0.0;
	double S = 0.0;
	double h, x;
	int flag;
	h = b - a;
	T1 = h / 2 * ( func( a ) + func( b ) );//计算第一个面积
	do {
		double s = 0;
		x = a + h /2;
		while( x < b )//循环计算
		{
			s = s + func( x );
			x = x + h;
		}
		T2 = T1 / 2 + h / 2 * s;
		if( fabs( T1 - T2 ) >= e )//判断精度
		{
			h = h / 2;
			T1 = T2;
			flag = 1;//当精度不够时，标志为1
		}
		else flag = 0;
	}
	while( flag );//为1时（意味着精度不够），继续计算。
	return T2;
}
//-------------------------------