为什么需要数值积分的C/C++实现
现在比较流行将Python作为科学计算语言来使用, 但是在使用Python的过程中有一些问题是无法或很难避开/解决的.
e.g. 程序运行速度稍慢(虽然说使用numpy可以改善); 由于存在CPython的全局解释器锁(GIL), 多线程反而比单线程要满; etc.
数值积分的C/C++实现
概述
首先我们先实现最基础的数值积分(暂时不考虑反常积分的问题)
数值积分的方法有很多, 这里笔者选择辛普森\frac{3}{8} 法则来计算
\int_{a}^{b} f(x)dx ≈ \frac{3h}{8} [f_0 + 3f_1 + 3f_2 + f_3], h = \frac{b-a}{3}
具体代码实现
long double integration(long double (*func)(long double), double a, double b, int n)
{
// Simpson’s Three-Eighths Rule
int count = 1;
float f, fact2 = 0, fact1 = 0, integral, fact3 = 0;
double h = (b - a) / n;
double i = a + h;
double f1 = func(a);
double f2 = func(a + n * h);
fact1 = f1 + f2;
if ((n % 3) != 0)
{
// Invalid No. of Subintervals!
return -1;
}
while (i <= b)
{
f = func(i);
if (count % 3 == 0)
{
if (i <= n - 3)
fact3 = fact3 + f;
}
else
{
if (i <= n - 1)
fact2 = fact2 + f;
}
i = i + h;
count++;
}
fact2 = fact2 * 3;
fact3 = fact3 * 2;
return (3 * h / 8) * (fact1 + fact2 + fact3);
}
请注意, 传入函数时无法多参数传入, 请有需求的读者自行修改一下
e.g.除积分变量x外另外传入一个参数p
long double integration_p2(long double (*func)(double, long double), double p, double a, double b, int n)
{
// Simpson’s Three-Eighths Rule
int count = 1;
float f, fact2 = 0, fact1 = 0, integral, fact3 = 0;
double h = (b - a) / n;
double i = a + h;
double f1 = func(p, a);
double f2 = func(p, a + n * h);
fact1 = f1 + f2;
if ((n % 3) != 0)
{
// Invalid No. of Subintervals!
return -1;
}
while (i <= b)
{
f = func(p, i);
if (count % 3 == 0)
{
if (i <= n - 3)
fact3 = fact3 + f;
}
else
{
if (i <= n - 1)
fact2 = fact2 + f;
}
i = i + h;
count++;
}
fact2 = fact2 * 3;
fact3 = fact3 * 2;
return (3 * h / 8) * (fact1 + fact2 + fact3);
}
读者比对可以发现, 其实就是在内部实现的时候向被积函数多传了一个参数
由于笔者刚从Python转到C/C++, 对其语法不太熟悉, 就不实现不定参数传入了(好多年前笔者明明学过C/C++的来着, 长时间不用都忘了... ...)
数值积分: 第一类反常积分
概述
无穷区间的积分又称第一类反常积分.常规计算方法是将积分上限+∞视为常数, 然后按照定积分来处理, 再将计算结果取极限.
计算第一类反常积分[图片上传失败...(image-5ffe6b-1660355781825)]
时, 可以计算[图片上传失败...(image-7a3c2f-1660355781825)]
和[图片上传失败...(image-e059b7-1660355781825)]
的差, 若两者差满足要求, 则停止运算. 若不满足要求, 则进一步迭代计算.
具体代码实现
unsigned long int upperLimit = 100;
long double IA = 0;
long double IB = 0;
while (abs(IA - IB) <= ESP || upperLimit / 100 >= CMAX)
{
IA = IB;
IB = integration(func, 0, upperLimit, N);
upperLimit = upperLimit + 100;
}
return 0;
这里的ESP是IA与IB的被允许最大差值, CMAX是最大迭代计算次数
数值积分: 第二类反常积分
概述
第二类反常积分是值积分区间包含奇异点(singular points). 常规计算方法是将积分积分区间在奇异点内收, 然后按照定积分来处理, 再将计算结果取极限.
这里直接引用了数值积分|第二类反常积分, 还请读者自行参考, 这里就不给出C/C++的实现了(有空再补上吧poi...)
算法
