通常大家都用微分法来求定积分的,但切割图形(微分)的方法也是有细微区别的。网上流行以下代码:
#include <iostream> #include <cmath> using namespace std; #define N 10000 #define Pi 3.1415926 double func(double x){ return sin(x); } double integral(double a, double b) { double r = 0, d = (b-a)/N; for(int i=0; i<N; i++) r += d * func(a+i*d); return r; } int main(void) { double result; result = integral(0, Pi/2); cout << "∫[0,Pi/2] SIN(x) = "; cout << result << endl; return 0; } /* ∫[0,Pi/2] SIN(x) = 0.999921 -------------------------------- Process exited after 0.828 seconds with return value 0 请按任意键继续. . . */
定积分的微分原理:
切割图形有四种方法,见上右图:红色矩形以左边为高,绿色矩形以右边为高,橙色矩形以中间值为高,第四种取红绿矩形的平均值即切割成梯形。
红色矩形在递增区间少算面积,递减则多算;绿色矩形反之。橙色矩形或梯形的误差稍微小些。下面用log(x)函数 [即以e为底的对数ln(),以10为底的函数名为log10()] 来比较一下,四种切割法的运算结果,为比较它们的误差把分割块数N值降至100:
#include <iostream> #include <cmath> using namespace std; #define N 100 #define E 2.7182818 double func(double x){ return log(x); } double integralL(double a, double b) { double r = 0, d = (b-a)/N; for(int i=0; i<N; i++) r += d * func(a+i*d); return r; } double integralR(double a, double b) { double r = 0, d = (b-a)/N; for(int i=0; i<N; i++) r += d * func(a+i*d+d); return r; } double integralM(double a, double b) { double r = 0, d = (b-a)/N; for(int i=0; i<N; i++) r += d * func(a+i*d+d/2); return r; } double integralT(double a, double b) { double r = 0, d = (b-a)/N; for(int i=0; i<N; i++) r += d * (func(a+i*d)+func(a+i*d+d))/2; return r; } int main(void) { double result; result = integralL(1, E); cout << "左高 ∫[1,E] LN(x)) = "; cout << result << endl; result = integralR(1, E); cout << "右高 ∫[1,E] LN(x) = "; cout << result << endl; result = integralM(1, E); cout << "中高 ∫[1,E] LN(x) = "; cout << result << endl; result = integralT(1, E); cout << "梯形 ∫[1,E] LN(x) = "; cout << result << endl; return 0; } /* 左高 ∫[1,E] LN(x) = 0.991393 右高 ∫[1,E] LN(x) = 1.00858 中高 ∫[1,E] LN(x) = 1.00001 梯形 ∫[1,E] LN(x) = 0.999984 -------------------------------- Process exited after 1.159 seconds with return value 0 请按任意键继续. . . */
从本例中看,高取中值的矩形在精度上明显胜出。以下切入正题:
函数作函数的形参
起因是要求很多个被积函数的定积分,如果不用函数作形参就要写很多个与每个被积函数一一对应的定积分函数。如:
#include <iostream> #include <cmath> using namespace std; #define N 100000 #define E 2.7182818 #define Pi 3.1415926 double func1(double x){ return sin(x); } double func2(double x){ return log(x); } double func3(double x){ return abs(x); } double func4(double x){ return exp(x); } double Integral1(double a, double b) { double r = 0, d = (b-a)/N; for(int i=0; i<N; i++) r+=d*func1(a+i*d+d/2); return r; } double Integral2(double a, double b) { double r = 0, d = (b-a)/N; for(int i=0; i<N; i++) r+=d*func2(a+i*d+d/2); return r; } double Integral3(double a, double b) { double r = 0, d = (b-a)/N; for(int i=0; i<N; i++) r+=d*func3(a+i*d+d/2); return r; } double Integral4(double a, double b) { double r = 0, d = (b-a)/N; for(int i=0; i<N; i++) r+=d*func4(a+i*d+d/2); return r; } int main(void) { double result=0; result = Integral1(0, Pi/2); cout << "∫[0, Pi/2] sin(x) =\t"; cout << result << endl; result = Integral2(1, E); cout << "∫[1, E] ln(x) =\t"; cout << result << endl; result = Integral3(-1, 1); cout << "∫[-1, 1] abs(x) =\t"; cout << result << endl; result = Integral4(-1, 1); cout << "∫[-1, 1] exp(x) =\t"; cout << result << endl; cout << "\n常数:e-1/e = " << E-1/E << endl; return 0; }
函数作参数就省事了,方法一:直接把被积函数声明作定积分函数的第三个参数
double Integral(double a, double b, double myfunc(double)) { double r = 0, d = (b-a)/N; for(int i=0; i<N; i++) r+=d*myfunc(a+i*d+d/2); return r; } // 调用时,只要用被积函数名替换声明的函数名myfunc即可: r1 = Integral(0, Pi/2, func1); r2 = Integral(1, E, func2); r3 = Integral(-1, 1, func3); r4 = Integral(-1, 1, func4);
方法二:定义一个函数指针,函数名即指针
double Integral(double a, double b, double (*fn)(double)) { double r = 0, d = (b-a)/N; for(int i=0; i<N; i++) r+=d*fn(a+i*d+d/2); // 注意:这里的*是乘号 return r; }
方法三:自定义一个函数数据类型
typedef double myFuncType(double); double Integral(double a, double b, myFuncType fn) { double r = 0, d = (b-a)/N; for(int i=0; i<N; i++) r+=d*fn(a+i*d+d/2); return r; }
方法四:自定义一个函数指针类型
typedef double (*FuncType)(double); double Integral(double a, double b, FuncType fn) { double r, d = (b-a)/N; for(int i=0; i<N; i++) r+=d*fn(a+i*d+d/2); return r; }
方法五:调用类的成员函数
class IntegralFunc { public: double fun1(double x){ return sin(x); } ... }; double CIntegral(double a, double b, double (IntegralFunc::*pf)(double)) { IntegralFunc fn; double r = 0, d = (b-a)/N; for(int i=0; i<N; i++) r+=d * (fn.*pf)(a+i*d+d/2); return r; } // 调用时: result = CIntegral(0, Pi/2, &IntegralFunc::fun1);
注意此方法中操作符和括号的用法,已着重加粗表示。
#include <iostream> #include <iomanip> #include <cmath> using namespace std; #define N 100000 #define E 2.7182818 #define Pi 3.1415926 class IntegralFunc { public: double fun1(double x){ return sin(x); } double fun2(double x){ return log(x); } double fun3(double x){ return abs(x); } double fun4(double x){ return exp(x); } }; double CIntegral(double a, double b, double (IntegralFunc::*pf)(double)) { IntegralFunc fn; double r = 0, d = (b-a)/N; for(int i=0; i<N; i++) r+=d * (fn.*pf)(a+i*d+d/2); return r; } int main(void) { double result=0; result = CIntegral(0,Pi/2,&IntegralFunc::fun1); cout << "∫[0, Pi/2] sin(x) =\t"; cout << result << endl; result = CIntegral(1, E, &IntegralFunc::fun2); cout << "∫[1, E] ln(x) =\t"; cout << result << endl; result = CIntegral(-1, 1, &IntegralFunc::fun3); cout << "∫[-1, 1] abs(x) =\t"; cout << result << endl; result = CIntegral(-1, 1, &IntegralFunc::fun4); cout << "∫[-1, 1] exp(x) =\t"; cout << result << endl; cout << "\n常数:e-1/e = " << E-1/E << endl; return 0; }
最后找些有点难度的定积分题来验证一下函数的精确度,见:《C++ 用自定义函数验证高等数学的定积分例题》。