标签：C++
#include<iostream>
#include<iomanip>
#include<cmath>
using namespace std;
//求积：梯形公式
double Func_Integral_Trapezoid
(double lo, double hi, double(*Func)(double), int n = 1000)
{//lo:下限；hi:上限；Func:函数；n:等分数
if (n <= 0) n = 100;
double x;
double step = (hi - lo) / n;
double result = 0.0;
x = lo;
for (int i = 1; i < n; i++) {
x += step;
result += Func(x);
}
result += (Func(lo) + Func(hi)) / 2;
result *= step;
return result;
}
//求积：Simpson公式
double Func_Integral_Simpson
(double lo, double hi, double(*Func)(double), int n = 1000)
{//lo:下限；hi:上限；Func:函数；n:等分数
if (n <= 0) n = 100;
double x;
double step = (hi - lo) / n;
double result1 = 0.0;
x = lo;
for (int i = 1; i < n; i++) {
x += step;
result1 += Func(x);
}
result1 *= 2;
double result2 = 0.0;
x = lo + step / 2;
for (int i = 0; i < n; i++) {
result2 += Func(x);
x += step;
}
result2 *= 4;
double result = result1 + result2 + Func(lo) + Func(hi);
result *= step / 6;
return result;
}
//求积：Cotes公式
double Func_Integral_Cotes
(double lo, double hi, double(*Func)(double), int n = 1000)
{//lo:下限；hi:上限；Func:函数；n:等分数
if (n <= 0) n = 100;
double x;
double step = (hi - lo) / n;
double result1 = 0.0;
x = lo;
for (int i = 1; i < n; i++) {
x += step;
result1 += Func(x);
}
result1 *= 14;
double result2 = 0.0;
x = lo + step / 2;
double result3 = 0.0;
double x1 = lo + step / 4;
double x2 = lo + step / 4 * 3;
for (int i = 0; i < n; i++) {
result2 += Func(x);
result3 += Func(x1) + Func(x2);
x += step;
x1 += step;
x2 += step;
}
result2 *= 12;
result3 *= 32;
double result4 = (Func(lo) + Func(hi)) * 7;
double result = result1 + result2 + result3 + result4;
result *= step / 90;
return result;
}
//求积：Romberg公式
double Func_Integral_Romberg
(double lo, double hi, double(*Func)(double), int k = 4)
{//lo:下限；hi:上限；Func:函数；k:等分指数(区间等分成2^k份)
int size = k + 1;
double *matrix = new double[size*size];
for (int i = 0; i < size*size; i++) matrix[i] = 0.0;
double step = hi - lo;
matrix[0] = Func_Integral_Trapezoid(lo, hi, Func, 1);
for (int i = 1; i < size; i++) {
int n = 1 << (i - 1);
for (int k = 0; k < n; k++) {
matrix[i*size + 0] += Func(lo + (k + 0.5)*step);
}
matrix[i*size + 0] *= step;
matrix[i*size + 0] += matrix[(i - 1)*size];
matrix[i*size + 0] /= 2.0;
step /= 2.0;
}
double temp = 1.0;
double factor1, factor2;
for (int j = 1; j < size; j++) {
temp *= 4.0;
factor1 = temp / (temp - 1);
factor2 = 1 / (temp - 1);
for (int i = j; i < size; i++) {
matrix[i*size + j] = factor1*matrix[i*size + j - 1]
- factor2*matrix[(i - 1)*size + j - 1];
}
}
double result = matrix[k*size + k];
delete[] matrix;
return result;
}
//测试用的被积函数，0到1积分为PI
double Func_test1(double x)
{
return 4 / (1 + x*x);
}
int main() {
cout << "Trapezoid Numerical Integration" << endl;
cout << setprecision(15) << Func_Integral_Trapezoid(0, 1, Func_test1, 1024) << endl << endl;
cout << "Simpson Numerical Integration" << endl;
cout << setprecision(15) << Func_Integral_Simpson(0, 1, Func_test1, 1024) << endl << endl;
cout << "Cotes Numerical Integration" << endl;
cout << setprecision(15) << Func_Integral_Cotes(0, 1, Func_test1, 1024) << endl << endl;
cout << "Romberg Numerical Integraion" << endl;
cout << setprecision(15) << Func_Integral_Romberg(0, 1, Func_test1, 10) << endl << endl;
getchar();
return 0;
} 标签：C++
原文地址：http://blog.51cto.com/gdutcxh/2116533
