数值计算：一重积分计算的C++实现

来源：互联网发布：法国制造业年度数据编辑：程序博客网时间：2024/05/22 19:26
#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积分为PIdouble 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;}
0 0