C++ 实现龙贝格算法

这里以一道例题为例

在这里为了方便计算，直接令

因为格式的原因，许多公式特殊符号复制不过来，成空格了。。

首先定义一个函数fx，方便后续的输入

#include<stdio.h>
#include<iostream>
#include<math.h>
#include <iomanip>
#define e        2.71828182845904523536
#define pi       3.14159265358979323846
using namespace std;

double fx(double  x) {		//令f(x)=	,求f(x)的值
	return 2 * pow(e, -1.0*x) / sqrt(pi);
}

然后，我的思路是先将 计算出，然后再据此算出

	t[0][0] = h/2.0*(fx(b)+fx(a));
	for (int i = 1; i < 10; i++)
	{
		double f = 0;
		x = a + h / 2.0;
		while (x < b) {
			f = f + fx(x);
			x += h;
		}
		t[0][i] = 0.5 * t[0][i - 1] + h/2 * f;  //递推公式求t(0)0
		h = h / 2.0;
	}

上面已经将t设置为全为0的10×10矩阵，现在只要判断t[i - 1][j - 1] 和 t[i - 1][j] 是否为0就可以知道能否求出t[i][j]的结果，书上的表是这样的：不过根据我的代码，上表的位置是t[1][1]而不是t[1][0]，以它所在的坐标来表示，代码如下：

for (int i = 1; i < 10; i++) {				//循环求
		for (int j = 1; j < 10; j++) {
			if (t[i - 1][j - 1] != 0 && t[i - 1][j] != 0) {	//判断能否求 的值
				t[i][j] = pow(4, i) / (pow(4, i) - 1) * t[i - 1][j] - t[i - 1][j - 1] / (pow(4, i) - 1);
				if (i == j) {
					if (abs(t[i][j] - t[i - 1][j - 1]) < pow(10, -7)) {	//检测是否达到精度，达到则停止循环并输出解
						cout << "卫星轨道的周长约为" << setprecision(8) << t[i][j]*4*7782.5<< "km" << endl;
						break;
					}
				}
			}
		}
	}

综合起来为：

#include<stdio.h>
#include<iostream>
#include<math.h>
#include <iomanip>
#define e        2.71828182845904523536
#define pi       3.14159265358979323846
using namespace std;

double fx(double  x) {		//令f(x)=	,求f(x)的值
	return 2 * pow(e, -1.0*x) / sqrt(pi);
}

int main()
{
	double t[10][10] = { 0 };			//创建一个二维数组
	double x; double h = 1;
	double a = 0, b = 1;		//a,b为上下限
	t[0][0] = h/2.0*(fx(1.0)+fx(0));		//求
	for (int i = 1; i < 10; i++)
	{
		double f = 0;
		x = a + h / 2.0;
		while (x < b) {
			f = f + fx(x);		//求
			x += h;
		}
		t[0][i] = 0.5 * t[0][i - 1] + h/2 * f;  //递推公式求
		h = h / 2.0;
	}
	for (int i = 1; i < 10; i++) {				//循环求
		for (int j = 1; j < 10; j++) {
			if (t[i - 1][j - 1] != 0 && t[i - 1][j] != 0) {	//判断能否求的值
				t[i][j] = pow(4, i) / (pow(4, i) - 1) * t[i - 1][j] - t[i - 1][j - 1] / (pow(4, i) - 1);
				if (i == j) {
					if (abs(t[i][j] - t[i - 1][j - 1]) < pow(10, -5)) {	
//检测是否达到精度，达到则停止循环并输出解
						cout << setprecision(8) << t[i][j] << "为所求解" << endl;
						break;
					}
				}
			}
		}
	}
}

代码其实是有些问题的，正常答案应该是类似这样的

但实际上这个代码已经从算到 了，不过最后的结果是一样的。