[numerical analysis exercise] numerical integration (complex trapezoid, complex Simpson, Romberg integral) C with STL implementation

Explain the algorithm a little

Complex trapezoid and complex Simpson In fact, it's all a train of thought , It's just that the compound trapezoid is added at one time 2 Two function values , Jump a space h distance ; And refolding Simpson It's a one-time addition 3 It's worth , Jump once 2h Distance of  .

The formula in the book is beautifully written , But it is not suitable for computer implementation , So the above idea is cyclic plus , And the book is one step .

Romberg The algorithm is troublesome , The formula given in the book is not clear enough , My code is an improvement on the formula in the book .

Logic is ：T Construction and SCR There are two kinds of constructions

T structure , The latter is equal to one-half of the former plus one delta_T

This delta_T yes , Sum of function values of all new points , Multiply by the new h It can be concluded that , The formula in the book is not suitable for computer implementation , So I modified it , It is also the first time in my life to use double Measure to write for loop .

SCR The structure belongs to the same system , Very simple formula .

The general idea is to divide the structure into two sections , First construct the first R1, So , We'll construct a total of T1-T8,S1-S4,C1-C2,R1, Such a triangular matrix .

Then there is a simple do while loop , Add accuracy control judgment .

Romberg Because it's an iterative method , Therefore, the accuracy can be approximately regarded as before and after R The difference between the

// Numerical integration 
// Complex trapezoidal integral and complex simpson The formula , Romberg formula 
#include<cstdio>
#include<climits>
#include<cmath>
#include<map> // Discrete storage 
double f(double x)
{
	return 1.0 / (x*x+1);
}
double trapzoid(double a,double b,int M)
{
	double ans = 0;
	double h = (b - a) / M;
	for (double x = a; x < b; x += h)
	{
		printf("%f x\n", x);
		ans += (f(x) + f(x + h));
	}
	ans = ans * h / 2; // The length of the middle segment of the compound trapezoid is h

	return ans;
}
double simpson(double a,double b,int M)
{
	if (M % 2)// Non even number 
	{
		puts(" Enter even number M");
		return -1;
	}
	
	double ans = 0;
	double h = (b - a) / M;
	for (double x = a; x < b; x += 2 * h)
	{
		printf("%f x\n", x);
		ans += (f(x) + 4 * f(x + h) + f(x + 2 * h));
	}
	ans = ans * 2 * h / 6; // The length of the compound Simpson section 2h

	return ans;
}
double romberg(double a, double b, double EPSILON)// Iteration is controlled by precision 
{
	std::map<int, double>T, S, C, R;
	T[1] = (b - a) * (f(a) + f(b)) / 2;
	printf("T[1] %f\n", T[1]);
	// Construct a triangular matrix 
	for (int Ti = 2; Ti <= 8; Ti *= 2)// Construct to T8
	{
		double delta_T = 0;
		double h = (b - a) / Ti;
		for (double start = a + h; start < b; start += 2 * h)// Insert the new node 
		{
			delta_T += f(start);
		}
		delta_T *= h;
		T[Ti] = T[Ti / 2] / 2 + delta_T;
		printf("T[%d] %f\n", Ti, T[Ti]);
	}
	for (int Si = 1; Si <= 4; Si *= 2)// structure Si To S4
	{
		S[Si] = pow(4,1) / (pow(4,1)-1) * T[Si * 2] 
			- 1.0 / (pow(4,1)-1) * T[Si];
		printf("S[%d] %f\n", Si, S[Si]);
	}
	for (int Ci = 1; Ci <= 2; Ci *= 2)// structure Ci To C2
	{
		C[Ci] = pow(4, 2) / (pow(4, 2) - 1) * S[Ci * 2] 
			- 1.0 / (pow(4, 2) - 1) * S[Ci];
		printf("C[%d] %f\n", Ci, C[Ci]);
	}
	R[1] = pow(4, 3) / (pow(4, 3) - 1) * C[2] // structure R1
		- 1.0 / (pow(4, 3) - 1) * C[1];
	printf("R[1] %f\n", R[1]);
	// Layer by layer structure R, End by precision control 
	int Ri = 1;
	do {
		Ri *= 2; 

		int Ti = Ri * 8;
		int Si = Ri * 4;
		int Ci = Ri * 2;
		double h = (b - a) / Ti;
		double delta_T = 0;
		// Current layer T
		for (double start = a + h; start < b; start += 2 * h)
		{
			delta_T += f(start);
		}
		delta_T *= h;
		T[Ti] = T[Ti / 2] / 2 + delta_T;
		printf("T[%d] %f\n", Ti, T[Ti]);
		// Current layer S
		S[Si] = pow(4, 1) / (pow(4, 1) - 1) * T[Si * 2]
			- 1.0 / (pow(4, 1) - 1) * T[Si];
		printf("S[%d] %f\n", Si, S[Si]);
		// Current layer C
		C[Ci] = pow(4, 2) / (pow(4, 2) - 1) * S[Ci * 2]
			- 1.0 / (pow(4, 2) - 1) * S[Ci];
		printf("C[%d] %f\n", Ci, C[Ci]);
		// Current layer R
		R[Ri] = pow(4, 3) / (pow(4, 3) - 1) * C[Ri * 2]
			- 1.0 / (pow(4, 3) - 1) * C[Ri];
		printf("R[%d] %f\n", Ri, R[Ri]);
	} while (R[Ri] - R[Ri / 2] > EPSILON);

	return R[Ri];
}

int main(void)
{
	double a, b, epsilon;
	int M;
	scanf("%lf %lf %d %lf", &a, &b, &M , &epsilon);
	printf("%f\n", trapzoid(a, b, M));
	printf("%f\n", simpson(a, b, M));
	printf("%f\n", romberg(a, b, epsilon));

	return 0;
}