数值计算中用复化梯形公式、复化抛物线公式及龙贝格求积公式求定积分的编程算法

三道c++编程题，分别用到复化梯形公式和复化辛普生公式、龙贝格算法、R-K算法对方程

第一题：

#include
#include
//复化梯形公式
doubleReiterationOfTrapezoidal(doublea,doubleb,doublen,doublef(doublex))
{
	doubleh,fa,fb,xk,t;
	h=(b-a)/n;//步长
	doubles=0.0;
	for(intk=1;k	{
		xk=a+k*h;
		s=s+f(xk);
	}
	returnh*(f(a)+f(b))/2+h*s;
}
//复化Simpson公式
doubleReiterationOfSimpson(doublea,doubleb,doublen,doublef(doublex))
{
	doubleh,fa,fb,xk,xj;
	h=(b-a)/n;//步长
	doubles1=0.0;
	doubles2=0.0;
	for(intk=1;k	{
		xk=a+k*h;
		s1=s1+f(xk);
	}
	for(intj=0;j	{
		xj=a+(j+0.5)*h;
		s2=s2+f(xj);
	}
	returnh/6*(f(a)+f(b)+2*s1+4*s2);
}
doublefunc(doublex)
{
	/*
	*1.0->10
	*1.1->11
	*...
	*2.0->20
	*/
	switch((int)round(x*10))
	{
	case10:return1;
	case11:return.90909;
	case12:return.83333;
	case13:return.76923;
	case14:return.71429;
	case15:return.66667;
	case16:return.625;
	case17:return.58824;
	case18:return.55556;
	case19:return.52632;
	case20:return.5;
	}
	return0;
}
intmain()
{
//都是分10步积完，从1到2
doubletra=ReiterationOfTrapezoidal(1,2,10,func);
doublesim=ReiterationOfSimpson(1,2,10,func);
printf("Trapezoidal:%f\n",tra);
printf("Simpson:%f\n",sim);
printf("Error:%f\n",(tra-sim));
return0;
}

第二题：

#include
#include
//龙贝格
doubleRomberg(doublea,doubleb,doubleeps,doublef(double))
{
	intk=1;
	doubleS,x,T1,T2,S1,S2,C1,C2,R1,R2,h=b-a;
	T1=h*(f(a)+f(b))/2;
	while(1)
	{
		S=0;x=a+h/2;
		do
		{
			S+=f(x);x+=h;
		}while(x		T2=(T1+h*S)/2.0;
		if(fabs(T2-T1)		S2=T2+(T2-T1)/3.0;
		if(k==1)
		{
			T1=T2;S1=S2;h/=2;k+=1;continue;
		}
		C2=S2+(S2-S1)/15.0;
		if(k==2)
		{
			C1=C2;T1=T2;S1=S2;h/=2.0;k+=1;continue;
		}
		R2=C2+(C2-C1)/63;
		if(k==3)
		{
			R1=R2;C1=C2;T1=T2;S1=S2;h/=2;k+=1;continue;
		}
		if(fabs(S1-S2)			return(S2);
		R1=R2;C1=C2;T1=T2;S1=S2;
		h/=2;
		k+=1;
		return(R1);
	}
}
doublefunc(doublex)
{
	doublec=cos(x);
	returnsqrt(x+c*c);
}
intmain()
{
	constdoublepi=3.14159265358979323;
	//从1到pi，精度是1e-5
	printf("%f\n",Romberg(0,pi,1e-5,func));
	return0;
}

MATLAB问题

不知道你的学历水平是什么样的，我尽量简单说

普通积分：

一般来说是找到被积函数的原函数，然后把上下限带入求值（定积分）

得到的解是精确解

参见大学课程：《高等数学》

在matlab中可用函数：int（），求解

数值积分：

由于在实际工程中绝大多数的积分都无法找到原函数

所以，想要得到精确解很困难

为了能够算出结果一些大牛（牛顿，高斯等等）搞出了数值解

即：数值积分（与精确解的误差满足需要）

具体的方法有很多比如：

梯形公式；辛普森公式；递归公式；龙贝格积分；自适应积分；高斯-勒让德积分等等

参见研究生课程：《数值分析》

在matlab中自带了一些求数值积分的函数：

trapz（）：基于复化梯形公式

integral（）：求解一元数值积分

integral2（）：求解一般区域二重积分数值解

integral3（）：求解一般区域三重积分数值解

我个人最常使用——高斯求积法求数值积分

一般来说取10-20个高斯点即可得到足够精度的数值解

最后：请叫我雷锋！

一道定积分？

题主给出的定积分问题，其被积函数比较复杂，所以用基本积分公式求解很难得到精确解。但我们可以通过数值积分的方法（如梯形公式，抛物线公式，龙贝格公式等等），得到其数值解。

对于题主的被积函数为x^x的定积分，使用复合抛物线公式，并借助于Excel来求解。

1、将积分限【1，2】分成10等份，即

k=0，1，2，3，。。。，10

2、计算积分计算点，即

x(k)=1，1.1，1.2，。。。，2

3、将被积函数写出 f(x)=x^x，并计算其y(k)值

4、根据复合抛物线公式计算，得到其积分值

5、计算过程如下

数值计算的构造数值积分

构造数值积分公式最通常的方法是用积分区间上的n 次插值多项式代替被积函数，由此导出的求积公式称为插值型求积公式。特别在节点分布等距的情形称为牛顿-柯茨公式，例如梯形公式与抛物线公式就是最基本的近似公式。但它们的精度较差。龙贝格算法是在区间逐次分半过程中，对梯形公式的近似值进行加权平均获得准确程度较高的积分近似值的一种方法，它具有公式简练、计算结果准确、使用方便、稳定性好等优点，因此在等距情形宜采用龙贝格求积公式。当用不等距节点进行计算时，常用高斯型求积公式计算，它在节点数目相同情况下，准确程度较高，稳定性好，而且还可以计算无穷积分。数值积分还是微分方程数值解法的重要依据。许多重要公式都可以用数值积分方程导出。

数值积分的基本思想

数学上，对于积分

地球物理数据处理基础

只要找到被积函数f（x）的原函数F（x），便可得到

地球物理数据处理基础

但在实际使用中，这种方法往往有困难，因为很多的被积函数找不到用初等函数表示的原函数，例如 等。在地球物理的实际问题中，被积函数f（x）往往不清楚或者很复杂，只是通过观测、测量等获取了一些离散的值yi＝f（x_i），即已知的f（x）是一个列表函数，无法求取其原函数。因此，这类问题就必须依靠数值积分解决。

对于列表函数f（x）的积分I^*＝∫^b_af（x）dx，当已知f（x）在区间［a，b］上n＋1个节点a＝x₀12<…n＝b的观测值f₀，f₁，f₂，…，f_n时，自然地想到用它们的插值多项式φ（x）的积分I＝∫^b_aφ（x）dx来近似代替I^*。其几何意义就是用［a，b］上插值曲线φ（x）与x轴所夹的面积代替f（x）曲线与x轴所夹的面积。

利用不同的插值多项式的积分可导出不同的求积系数和求积公式，通常利用分段低阶的插值多项式求积应用最广，表7－1给出了不同数值积分方法的优缺点。复化辛卜生公式和复化梯形公式既简便易行、容易理解，又具有较高的精度，因此是地球物理计算中的主要的求积方法。高斯公式选点少精度高，是正演计算的求积方法。至于龙贝格算法，只要能加密节点，就可以达到最快的收敛速度，且每次提高二阶逼近的精度，是数值积分中很巧妙、很优秀的方法，但是由于地球物理测量的成本和周期限制观测点密度是有限的，因而无法发挥这种算法的优势。样条求积精度高，是地球物理计算中理想的求积方法。基于这些原因本章重点介绍复化梯形求积、复化辛卜生求积，至于样条求积，在第六章样条插值函数的基础上很容易实现。

表7－1 不同求积公式对比表