数值积分: 梯形规则--复合梯形规则--辛普森规则--复合辛普森规则--龙贝格求积公式

本文介绍了数值积分中的几种重要方法,包括梯形规则、辛普森规则及其复合形式,以及龙贝格求积公式。通过这些方法可以有效解决复杂的积分问题。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

数值积分:梯形规则--复合梯形规则--辛普森规则--复合辛普森规则--龙贝格求积公式


1.问题描述

微积分方法求积有很大的局限性,当碰到被积函数很复杂时,找不到相应的原函数。积分值

在几何上可解释为由 x=a,x=b,y=0和y=f(x) 所围成的曲边梯形的面积。积分计算之所以有困难,就是因为这个曲边梯形有一条边y=f(x)是曲线。


2.理论与方法

依据积分中值定理,底为b-a,而高为f(e)的矩形面积恰等于所求曲边梯形的面积I.f(e)称作区间[a,b]上的平均高度。这样,只要对平均高度f(e)提供一种算法,便相应地获得一种数值求积的算法。

1.梯形规则(Trapezoidal rule)
简单选取区间[a ,b]的中点高度作为平均高度。取h=b-a
        a0=(a-b)(x-b)/(a-b)dx=(b-a)/2
        a1=(a-b)(x-a)/(b-a)dx=(b-a)/2
得到:
2.辛普森规则(Simpson rule)

可视作用a , bc=(a+b)/2三点高度的加权平均值作为平均高度。

3.复合梯形规则(Composite numerical)
设将求积区间[a,b]划分为n等份,步长h=(b-a)/2 ,等分点为xi=a+bi , i=0,1,...,n 所谓复化求积法,就是先用低阶求积公式求得每个子段[xi,xi+1]上的积分值,然后再将它们累加求和,用各段积分之和Ii,i=0,1,n-1作为所求积分的近似值。
复化梯形公式:
4.复合辛普森规则(Composite Simpson)

记子段[xi,xi+1]的中点为则复化公式为
复化Simpson公式

5.龙贝格求积公式(Romberg)
龙贝格求积公式也称为逐次分半加速法。它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上,构造出一种加速计算积分的方法。 作为一种外推算法, 它在不增加计算量的前提下提高了误差的精度.

在等距基点的情况下,用计算机计算积分值通常都采用把 区间 逐次分半的方法进行。这样,前一次分割得到的函数值在分半以后仍可被利用,且易于编程。

对区间[a b],令h=b-a构造梯形值序列{T2k}

T1=h[f(a)+f(b)]/2

把区间二等分,每个小区间长度为 h/2=(b-a)/2,于是

T2 =/2+[h/2]f(a+h/2

把区间四(2)等分,每个小区间长度为h/4 =b-a/4,于是

T4 =/2+[h/2][fa+h/4+fa+3h/4).....................
[ab] 2等分,分点=a+b-a/ 2 ·i i =012 · · · 2k)每个小区间长度为(b-a/ 2 .

3.算法与设计
1. 梯形规则 ( Trapezoidal rule )
计算区间端点函数值,代入公式:
2. 辛普森规则( Simpson rule )
计算区间端点及中点函数值,代入公式

3.复合梯形,辛普森规则(Composite)

4.龙贝格求积公式(Romberg)
已给函数f(x),求积区间[a,b] 及精度e

准备初值:令 h=b-a,n=1 在步长二分过程中,逐步生成,并存数据T8,S4,C2,R1于单元t1,s1,c1,r1
二分求梯形值 h=h/2,n=2n 根据公式求出新的梯形值,存于单元t2

松弛加速:根据公式求加速值,存于单元s2,c2,r2

精度控制:如果|r2-r1|<e,则输出r2作为结果,终止计算;
否则t1=t2,s1=s2,c1=c2,r1=r2,转2继续二分。

4.代码
void tra()
{
	cout << endl;
	cout << "trapezoidal rule" << endl;
	double a = 0, b = 0.8;
	double fa = 0.2 + 25 * a - 200 * pow(a,2) + 675 * pow(a, 3) - 900 * pow(a, 4) + 400 * pow(a, 5);
	double fb = 0.2 + 25 * b - 200 * pow(b, 2) + 675 * pow(b, 3) - 900 * pow(b, 4) + 400 * pow(b, 5);
	double n = 100;
	double h = (b - a) / n;
	double w[100]; w[0] = 0;
	double f[100];
	for (int i = 1; i < n; i++) w[i] = w[i - 1] + h;
	for (int i = 1; i < n; i++)
	{
		f[i] = 0.2 + 25 * w[i] - 200 * pow(w[i], 2) + 675 * pow(w[i], 3) - 900 * pow(w[i], 4) + 400 * pow(w[i], 5);
	}
	double s = 0;
	s = ((h - a) / 2)*(fa + f[1]);
	for (int i = 2; i < n; i++)
	{
		s += ((h / 2)*(f[i] + f[i - 1]));
	}
	s+= ((h / 2)*(fb + f[99]));
	cout << s << endl;
}

void sim()
{
	cout << endl;
	cout << "simpson's rule" << endl;
	double a = 0, b = 0.8;
	double fa = 0.2 + 25 * a - 200 * pow(a, 2) + 675 * pow(a, 3) - 900 * pow(a, 4) + 400 * pow(a, 5);
	double fb = 0.2 + 25 * b - 200 * pow(b, 2) + 675 * pow(b, 3) - 900 * pow(b, 4) + 400 * pow(b, 5);

	double n = 100;
	double h = (b - a) / n;
	double w[100]; w[0] = 0;
	double f[100];
	for (int i = 1; i < n; i++) w[i] = w[i - 1] + h;
	for (int i = 1; i < n; i++)
	{
		f[i] = 0.2 + 25 * w[i] - 200 * pow(w[i], 2) + 675 * pow(w[i], 3) - 900 * pow(w[i], 4) + 400 * pow(w[i], 5);
	}
	double s = 0;
	f[0] = fa;

	for (int i = 1; i < n-1; i+=2)
	{
		s += (2*h / 6)*(f[i-1] + 4*f[i]+f[i+1]);
	}

	cout << s << endl;
}

void com_tra()
{
	cout << endl;
	cout << "the composite trapezoidal rule" << endl;
//	double h = 0.1;
	double a = 0, b = 0.8;
	double fa = 0.2 + 25 * a - 200 * pow(a, 2) + 675 * pow(a, 3) - 900 * pow(a, 4) + 400 * pow(a, 5);
//	double f7 = 0.2 + 25 * h*7 - 200 * pow(h * 7, 2) + 675 * pow(h * 7, 3) - 900 * pow(h * 7, 4) + 400 * pow(h * 7, 5);
	double fb = 0.2 + 25 * b - 200 * pow(b, 2) + 675 * pow(b, 3) - 900 * pow(b, 4) + 400 * pow(b, 5);
	
	double n = 100;
	double h = (b - a) / n;
	double w[100]; w[0] = 0;
	double f[100];
	for (int i = 1; i < n; i++) w[i] = w[i - 1] + h;
	for (int i = 1; i < n; i++)
	{
		f[i] = 0.2 + 25 * w[i] - 200 * pow(w[i], 2) + 675 * pow(w[i], 3) - 900 * pow(w[i], 4) + 400 * pow(w[i], 5);
	}
	double s = 0;
	f[0] = fa;
	for (int i = 1; i < n ; i++)
	{
		s += (h / 2)*( 2 * f[i]);
	}
	s += (h / 2)*(fa + fb);
	cout << s << endl;
}

void com_sim()
{
	cout << endl;
	cout << "the composite simpson's rule" << endl;
	//	double h = 0.1;
	double a = 0, b = 0.8;
	double fa = 0.2 + 25 * a - 200 * pow(a, 2) + 675 * pow(a, 3) - 900 * pow(a, 4) + 400 * pow(a, 5);
	double fb = 0.2 + 25 * b - 200 * pow(b, 2) + 675 * pow(b, 3) - 900 * pow(b, 4) + 400 * pow(b, 5);

	double n = 100;
	double h = (b - a) / n;
	double w[100]; w[0] = 0;
	double f[100];
	for (int i = 1; i < n; i++) w[i] = w[i - 1] + h;
	for (int i = 1; i < n; i++)
	{
		f[i] = 0.2 + 25 * w[i] - 200 * pow(w[i], 2) + 675 * pow(w[i], 3) - 900 * pow(w[i], 4) + 400 * pow(w[i], 5);
	}
	double s = 0;
	f[0] = fa;
	for (int i = 1; i < n; i++)
	{
		if (i % 2) s += (h / 3)*(2 * f[i]);
		else s += (h / 3)*(4 * f[i]);
	}
	s += (h / 3)*(fa + fb);
	cout << s << endl;
}


End
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值