龙贝格算法的实现以及与复合梯形公式精度的比较

一、实验名称：
数值积分的实现
二、实验目的：
a.验证龙贝格公式的积分效果
b.验证龙贝格公式对复合梯形公式精度的提高
三、实验内容
（一）用龙贝格公式求定积分
1.用户输入积分区间端点a，b，输入所求积分的精度值precision。
2.求初始步长h=b-a，此时的步长为第0次二分的结果。
3.求二分0次时的梯形公式积分结果T1.
4.进行第一次二分，求此时的梯形积分结果T2。
5.由已经求出的T1和T2计算出辛普生积分结果的S2.
6.继续计算梯形积分公式的积分结果，以此递推求龙贝格积分结果。
7.判断此时的龙贝格积分结果精度是否达到要求，若未达到要求，则继续二分计算，否则输出积分结果。
（二）验证不同积分公式误差的大小
1.用户输入积分区间端点a，b，输入梯形公式二分的次数n。
2.定义数组T，S，C，R分别用来存储梯形，辛普生型，科特斯型，龙贝格型积分结果。
3.调用梯形积分函数，计算T。
4.调用龙贝格积分函数，进行递归计算R。
5.分别输出上述四种类型的积分，同时，每种积分类型的输出中都包括积分结果，误差，有效数字，与上次误差与本次误差的比值。
四、程序关键语句描述
（一）用龙贝格公式求定积分

double f(double x) {
	if (x == 0)
		return 1;
	else return sin(x) / x;
}

上述代码是用来计算函数值f（x），这样就不用每次都把节点值输入到计算机了。判断x是否等于0的原因是计算机不会求极限，因此需要考虑特殊点的函数值。

	h = b - a;
	T1 = h / 2.0 * (f(a) + f(b));
	k = 1;

计算出h的初始步长，T1为二分0次的梯形复合梯形积分结果，k用来判断进行了几次二分，因为至少要二分三次才能计算龙贝格积分结果。

while (1) {
		S = 0;
		x = a + h / 2.0;
		while (1) {	//计算从i = 1到i = n-1的f（xi）的和，并且用S存储
			S = S + f(x);
			x = x + h;
			if (x >= b)	//若已超过区间右端点
				break;
		}
		T2 = T1 / 2.0 + h / 2.0 * S;
		S2 = T2 + 1.0 / 3.0 * (T2 - T1);
		if (k == 1) {
			k++;
			h = h / 2.0;
			T1 = T2;
			S1 = S2;
			continue;
		}
		C2 = S2 + 1.0 / 15.0 * (S2 - S1);
		if (k == 2) {
			k++;
			h = h / 2.0;
			C1 = C2;
			T1 = T2;
			S1 = S2;
			continue;
		}
		R2 = C2 + 1.0 / 63.0 * (C2 - C1);
		if (k == 3) {
			k++;
			h = h / 2.0;
			R1 = R2;
			C1 = C2;
			T1 = T2;
			S1 = S2;
			continue;
		}
		if (fabs(R2 - R1) >= precision) {
			k++;
			h = h / 2.0;
			R1 = R2;
			C1 = C2;
			T1 = T2;
			S1 = S2;
			continue;
		}
		break;
	}

上面的代码中，共有两个while循环，内层while循环用来计算从i = 1到i = n-1的f（xi）的和，并且用S存储，当k=1时，让h减半，也就相当于进行了一次二分，同样的道理，k=2和k=3时进行相同的操作，当k=4时，现在已经计算出来了R2和R1的值，两者做差，判断精度是否满足要求，若不满足则继续二分下去，否则break语句跳出循环。
（三）验证不同积分公式误差的大小

double f(double x) {
	return pow(x,4) + pow(x,3) + pow(x,2) + pow(x,1) + 3 * pow(x,5) + 6 * pow(x,6);
}

上述代码用来计算函数值，本次实验所选取的函数为
f(x) = x + x2 + x3 + x4 + 3x5 + 6x6。
这个函数具有足够的复杂度，而且能够求出原函数，能够明显的观察龙贝格算法对梯形复合公式精度的提高。本实验选取的积分区间为[0,1],积分结果的真实值为1109/420。

double sum(int n, double h) {
	double xi;
	double result = 0;
	for (int i = 1; i <= n - 1; i++) {
		xi = a + i * h;
		result  = result + f(xi);
	}
	return result;
}

这是一个求和函数，在用复合梯形公式求积分的时候会用到，返回的结果为从第一个函数值到第n-1个函数值的累加和。

int valid_num(double r) {
	int l = 10;
	while (1) {
		if (r <= 0.5 * pow(10, 1-l) || l == 0)
			return l;
			l--;
	}
}

计算有效数字的函数，如果积分结果等于真实值，这时的有效数字从理论上来讲应该是正无穷，但显然计算机无法计算正无穷，因此假定有效数字l最高为十位。

void trapezoid() {	//梯形复合公式
	double h;
	for (int i = 0; i < k + 1; i++) {
		h = (b - a) / pow(2, i);
		T[i] = h / 2 * (f(a) + f(b) + 2 * sum(pow(2,i),h));
	}
}

计算复合梯形公式的函数，for循环中，i表示二分的次数，因此pow(2,i)相当于二分出来的小区间的个数n。

void Romberg() {	//龙贝格公式
	int i;
	for (i = 0; i < k; i++) 
		S[i] = 4.0 / 3.0 * T[i + 1] - 1.0 / 3.0 * T[i];
	for (i = 0; i < k - 1; i++)
		C[i] = 16.0 / 15.0 * S[i + 1] - 1.0 / 15.0 * S[i];
	for (i = 0; i < k - 2; i++)
		R[i] = 64.0 / 63.0 * C[i + 1] - 1.0 / 63.0 * C[i];
}

用龙贝格积分公式计算积分值的函数，公式来源于课本。

	for (i = 0,flag = 0; i < k + 1; i++,flag++) {
		printf("现在输出第%d次梯形积分\n",i);
		printf("积分值为%.20lf\n",T[i]);
		printf("误差为%.20lf\n", fabs(T[i] - 1109.0 / 420.0));
		printf("有效数字为%d\n", valid_num(fabs(T[i] - 1109.0 / 420.0)));
		if (flag != 0)
			printf("上次误差与本次误差的比值%lf\n",ter / fabs(T[i] - 1109 / 420.0));
		printf("\n");
		ter = fabs(T[i] - 1109.0 / 420.0);
	}

以梯形积分公式的输出为例，来说明一下输出过程的实现，double型变量ter保存的是本次积分结果的误差，目的是能够方便的计算相邻两次积分误差的比值，便于观察二分次数以及不同积分公式的误差情况。

六、实验体会
通过本次实验，用C语言实现了龙贝格积分公式，成功的验证了龙贝格积分公式的正确性，同时认识到此公式的巧妙之处。在编写程序的时候，涉及到的循环跳转比较多，刚开始感觉很混乱，后来画了一下流程图，感觉豁然开朗，用continue和break进行循环的控制，简化了代码量，同时增强了程序的可读性。
在验证龙贝格公式对梯形复合公式的提高时，我同时验证了几种不同积分公式的余项定理，梯形公式的余项是与h的平方成正比，因此相邻两次梯形积分的误差的比值应接近于4，同理，辛普生余项与h的四次方成正比，相邻两次辛普生积分的误差的比值应接近16，科特斯余项与h的六次方成正比，相邻两次科特斯积分的误差的比值应接近64，通过实验验证，上述理论是正确的。当进行三次二分时，从输出结果中能够看到，此时的龙贝格积分结果的精度比进行十次二分时的梯形复合积分结果的精度还要高，这从侧面说明了，龙贝格公式可以用尽量少的计算，得到较高的精度。同时，我也感受到了龙贝格公式的智慧所在。当进行十次二分时，我们发现第七次科特斯积分结果的输出中，误差比值变为37，第八次科特斯积分结果的输出中，误差比值变为0.6666，前面其余的结果都是符合理论的，只有这两个结果偏离理论，初步分析，应该是随着二分次数的增多，误差呈指数的趋势减小，而计算机又不能进行无限运算，因此怀疑是计算机本身的精度造成了这一现象。这也提示我们，在实际应用中，不能只考虑理论结果，还要考虑运算工具的运算能力，这样才能得到最佳的效率比。
最后，附上代码

//用龙贝格公式验证积分结果的正确性
#include "pch.h"
#include<stdio.h>
#include<math.h>
double f(double x) {
	if (x == 0)
		return 1;
	else return sin(x) / x;
}
int main()
{
	double a, b, precision, h, x;//a,b为区间的端点，presicion为给定的积分精度，h是每个积分小区间的长度
	double S;
	double T1, T2, S1, S2, C1, C2, R1, R2;
	int k;	//
	printf("请输入区间端点a，b\n");
	scanf("%lf%lf", &a, &b);
	printf("请输入精度\n");
	scanf("%lf", &precision);
	h = b - a;
	T1 = h / 2.0 * (f(a) + f(b));
	k = 1;
	while (1) {
		S = 0;	
		x = a + h / 2.0;
		while (1) {	//计算从i = 1到i = n-1的f（xi）的和，并且用S存储
			S = S + f(x);
			x = x + h;
			if (x >= b)	//若已超过区间右端点
				break;
		}
		T2 = T1 / 2.0 + h / 2.0 * S;
		S2 = T2 + 1.0 / 3.0 * (T2 - T1);
		if (k == 1) {
			k++;
			h = h / 2.0;
			T1 = T2;
			S1 = S2;
			continue;
		}
		C2 = S2 + 1.0 / 15.0 * (S2 - S1);
		if (k == 2) {
			k++;
			h = h / 2.0;
			C1 = C2;
			T1 = T2;
			S1 = S2;
			continue;
		}
		R2 = C2 + 1.0 / 63.0 * (C2 - C1);
		if (k == 3) {
			k++;
			h = h / 2.0;
			R1 = R2;
			C1 = C2;
			T1 = T2;
			S1 = S2;
			continue;
		}
		if (fabs(R2 - R1) >= precision) {
			k++;
			h = h / 2.0;
			R1 = R2;
			C1 = C2;
			T1 = T2;
			S1 = S2;
			continue;
		}
		break;
	}
	printf("精度为%lf的积分结果为：%.9lf\n", precision, R2);
}

 // 比较龙贝格公式对积分精度的提高效果
#include "pch.h"
#include<stdio.h>
#include<math.h>

int k;
double a, b;
double T[11];
double S[11];
double C[10];
double R[10];
double f(double x) {
	return pow(x,4) + pow(x,3) + pow(x,2) + pow(x,1) + 3 * pow(x,5) + 6 * pow(x,6);
}
double sum(int n, double h) {
	double xi;
	double result = 0;
	for (int i = 1; i <= n - 1; i++) {
		xi = a + i * h;
		result  = result + f(xi);
	}
	return result;
}
int valid_num(double r) {
	int l = 10;
	while (1) {
		if (r <= 0.5 * pow(10, 1-l) || l == 0)
			return l;
			l--;
	}
}
void trapezoid() {	//梯形复合公式
	double h;
	for (int i = 0; i < k + 1; i++) {
		h = (b - a) / pow(2, i);
		T[i] = h / 2 * (f(a) + f(b) + 2 * sum(pow(2,i),h));
	}
}
void Romberg() {	//龙贝格公式
	int i;
	for (i = 0; i < k; i++) 
		S[i] = 4.0 / 3.0 * T[i + 1] - 1.0 / 3.0 * T[i];
	for (i = 0; i < k - 1; i++)
		C[i] = 16.0 / 15.0 * S[i + 1] - 1.0 / 15.0 * S[i];
	for (i = 0; i < k - 2; i++)
		R[i] = 64.0 / 63.0 * C[i + 1] - 1.0 / 63.0 * C[i];
}
int main()
{
	printf("请输入a，b\n");
	scanf("%lf%lf", &a, &b);
	printf("请输入二分次数k(k > 2)\n");
	scanf("%d", &k);
	trapezoid();
	Romberg();
	int i,flag;
	double ter,rate;
	printf("\n");
	//printf("T\n");
	for (i = 0,flag = 0; i < k + 1; i++,flag++) {
		printf("现在输出第%d次梯形积分\n",i);
		printf("积分值为%.20lf\n",T[i]);
		printf("误差为%.20lf\n", fabs(T[i] - 1109.0 / 420.0));
		printf("有效数字为%d\n", valid_num(fabs(T[i] - 1109.0 / 420.0)));
		if (flag != 0)
			printf("上次误差与本次误差的比值%lf\n",ter / fabs(T[i] - 1109 / 420.0));
		printf("\n");
		ter = fabs(T[i] - 1109.0 / 420.0);
	}
	//printf("S\n");
	for (i = 0,flag = 0; i < k; i++,flag++) {
		printf("现在输出第%d次辛普生积分\n", i);
		printf("积分值为%.20lf\n", S[i]);
		printf("误差为%.20lf\n", fabs(S[i] - 1109 /420.0)); 
		printf("有效数字为%d\n", valid_num(fabs(S[i] - 1109.0 / 420.0)));
		if (flag != 0)
			printf("上次误差与本次误差的比值%lf\n", ter / fabs(S[i] - 1109 / 420.0));
		printf("\n");
		ter = fabs(S[i] - 1109.0 / 420.0);
	}
	//printf("C\n");
	for (i = 0,flag = 0; i < k - 1; i++,flag++) {
		printf("现在输出第%d次科特斯积分\n", i);
		printf("积分值为%.20lf\n", C[i]);
		printf("误差为%.20lf\n", fabs(C[i] - 1109 / 420.0));
		printf("有效数字为%d\n", valid_num(fabs(C[i] - 1109.0 / 420.0)));
		if (flag != 0)
			printf("上次误差与本次误差的比值%lf\n", ter / fabs(C[i] - 1109 / 420.0));
		printf("\n");
		ter = fabs(C[i] - 1109.0 / 420.0);
	}
	//printf("R\n");
	for (i = 0; i < k - 2; i++) {
		printf("现在输出第%d次龙贝格积分\n", i);
		printf("积分值为%.10lf\n", R[i]);
		printf("误差为%.10lf\n", fabs(R[i] - 1109 / 420.0));
		printf("有效数字为%d\n", valid_num(fabs(R[i] - 1109.0 / 420.0)));
	}
}

其实，用龙贝格公式验证积分结果的正确性还有一种写法，可以用goto语句实现，代码如下。这个代码可以对照着文末的流程图来看，这样会对这个程序理解的更深入。

// 龙贝格算法2.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//

#include "pch.h"
#include<stdio.h>
#include<math.h>
/*
int k;
double a, b;
double T[11];
double S[11];
double C[10];
double R[10];
double f(double x) {
	if (x == 0)
		return 1;
	return sin(x) / x;
}
double sum(int n, double h) {
	double xi;
	double result = 0;
	for (int i = 1; i <= n - 1; i++) {
		xi = a + i * h;
		result = result + f(xi);
	}
	return result;
}
void trapezoid() {	//梯形复合公式
	double h;
	for (int i = 0; i < k + 1; i++) {
		h = (b - a) / pow(2, i);
		T[i] = h / 2 * (f(a) + f(b) + 2 * sum(pow(2, i), h));
	}
}
void Romberg() {	//龙贝格公式
	int i;
	for (i = 0; i < k; i++)
		S[i] = 4.0 / 3.0 * T[i + 1] - 1.0 / 3.0 * T[i];
	for (i = 0; i < k - 1; i++)
		C[i] = 16.0 / 15.0 * S[i + 1] - 1.0 / 15.0 * S[i];
	for (i = 0; i < k - 2; i++)
		R[i] = 64.0 / 63.0 * C[i + 1] - 1.0 / 63.0 * C[i];
}
*/
double f(double x) {
	if (x == 0)
		return 1;
	return sin(x) / x;
}
int main()
{
	double a, b,r,h,x;
	double S;
	double T1, T2, S1, S2,C1,C2,R1,R2;
	int k;
	printf("请输入区间端点a，b\n");
	scanf("%lf%lf", &a, &b);
	printf("请输入精度\n");
	scanf("%lf", &r);
	h = b - a;
	T1 = h / 2.0 * (f(a) + f(b));
	k = 1;
L1:
	S = 0;
	x = a + h / 2.0;
L2:
	S = S + f(x);
	x = x + h;
	if (x < b)
		goto L2;
	T2 = T1 / 2.0 + h / 2.0 * S;

	S2 = T2 + 1.0 / 3.0 * (T2 - T1);
	if (k == 1) {
		k++;
		h = h / 2.0;
		T1 = T2;
		S1 = S2;
		goto L1;
	}
	C2 = S2 + 1.0 / 15.0 * (S2 - S1);
	if (k == 2) {
		C1 = C2;
		k++;
		h = h / 2.0;
		T1 = T2;
		S1 = S2;
		goto L1;
	}
	R2 = C2 + 1.0 / 63.0 * (C2 - C1);
	if (k == 3) {
		R1 = R2;
		C1 = C2;
		k++;
		h = h / 2.0;
		T1 = T2;
		S1 = S2;
		goto L1;
	}
	if (fabs(R2 - R1) >= r) {
		R1 = R2;
		C1 = C2;
		k++;
		h = h / 2.0;
		T1 = T2;
		S1 = S2;
		goto L1;
	}
	printf("精度为%lf的积分结果为：%.7lf\n", r,R2);
}