C语言积分思想和算法学习

C语言积分思想

标准C语言的积分好像很难，至少收集资料前是这么想的，当收集完资料，看到各位前辈的介绍，又重拾信心，梳理记录一下，以备后用吧。翻了翻参考文献，打底利用数值计算的方法都是利用微积分的思维，小图像累加求和来做的。实用的算法分为梯形算法（Trapezoidal Rule）和辛普森算法（Simpson’s Rule），那么区别和性能呢？本文的图形都来自文献Integration

Trapezoidal Rule

相对于上图的样条矩形，肉眼可见的误差很大，真实面积更大，如何减小误差，直接连接插值点，矩形就变成梯形了
很明显，误差小了，计算公式$$S=\frac{\delta x}{2}(y_0+y_1+y_1+y_2...+y_4+y_5)$$
引申一下，假设这个函数为$f(x)$积分范围$a-b$，把这个积分面积分成$n$份，采样点$n+1$。则积分（面积）可以表达为：$${\int }_a^{b}f(x)=\frac{b-a}{2\ast n}(f(a)+2\ast \sum _{k=1}^{n-1}f(x_k)+f(b))$$

Simpson’s Rule

这个规则就是认为利用抛物线来拟合边界会带来更好的结果，

$${\int }_a^{b}f(x)=\frac{\delta x}{3}(y_0+4y_1+2y_2+4y_3+2y_4+4y_5+y_6)$$
这个规则要求$n$必须是偶数，一般公式为：$${\int }_a^{b}f(x)=\frac{\delta x}{3}(y_0+2\sum y_{even}+4\sum y_{old}+y_n)$$
看到这个公式有点懵，不知从哪里来的，那先了解一下抛物线（parabolas）的面积公式，有幸搜到一篇博文可视化微积分：一种全新的视野解释抛物线的面积，大致上就是说围成抛物线的面积等于外围矩形面积的$1/3or2/3$，那要看怎么看了，上面就是把每条分成矩形和一个抛物线形，两个面积累加，上式采用$1/3stepto2/3$ 的方法，得到了这个规则公式，推导上图片吧，不码公式了。
这就提供了简单的积分运算的数值方法。可以利用c最基本的库来实现积分了。当然还是有误差的，参考文档里有几个现成的实例。

实现中等分问题的取舍：变步长和预制误差

前面的章节说明了数值积分的基本思想，但是应用中首先会遇到两个问题：

1. 等分多少合适？
2. 计算误差能符合要求吗？

其实这两个问题是有联系的，等分越多，误差越小。基于此，如果事先定好能接受的误差，用迭代的方法，发现两次等分结果的误差已经收敛到这个范围，那么一般认为结果已经满足要求了，关于这个算法背后的数学推导，暂时不去深究，首先了解一下变步长梯形算法。

变步长梯度求积

变步长即根据预定额误差决定等分多少，到达误差要求返回结果，假设待积分函数为$$T={\int }_a^{b}f(x)dx$$求解步骤：

1. 首先预定一个收敛误差$ϵ$步长为1：求解$f(a)$和$f(b)$围成的梯形，梯形的高为$b-a$,即$$\begin{gathered} T_1=\frac{b-a}{2}(f(a)+f(b)) \\ =\frac{1}{2}(\frac{b-a}{2}(f(a)+f(b)))+\frac{b-a}{2}(f(a)+f(b)) \end{gathered}$$
2. 然后步长为2：$$\begin{gathered} T_2=\frac{b-a}{4}(f(a)+2\ast f(a+b)+f(b)) \\ =\frac{1}{2}{T}_1+\frac{b-a}{2}f(a+b) \end{gathered}$$两次积分结果的误差$$|T_1-T_2|<ϵ$$，则完美的达到要求，可以返回$T_2$，否则还得继续迭代下去。
3. 迭代之前如果能找到一些减小运算量的办法最好，回顾梯形的一般公式： T n = b − a 2 ∑ k = 0 n − 1 { f ( x k ) + f ( x k + 1 ) } , x 0 = a , x n = b T_n=\frac{b-a}{2}\sum_{k=0}^{n-1}\{f(x_k)+f(x_{k+1})\}, x_{0}=a, x_n=b Tn​=2b−a​k=0∑n−1​{f(xk​)+f(xk+1​)},x0​=a,xn​=b那么如果继续二等分的话， T 2 n = b − a 2 ∑ k = 0 n − 1 { f ( x k ) + f ( ( x k + x k + 1 ) / 2 ) 2 + f ( ( x k + x k + 1 ) / 2 ) + f ( x k + 1 ) 2 } = b − a 4 ∑ k = 0 n − 1 { f ( x k ) + f ( x k + 1 ) } + b − a 2 ∑ k = 0 n − 1 f ( ( x k + x k + 1 ) / 2 ) = 1 2 T n + b − a 2 ∑ k = 0 n − 1 f ( ( x k + x k + 1 ) / 2 ) T_{2n}=\frac{b-a}{2}\sum_{k=0}^{n-1}\{\frac{f(x_k)+f((x_k+x_{k+1})/2)}{2}+\frac{f((x_k+x_{k+1})/2)+f(x_{k+1})}{2}\}\\=\frac{b-a}{4}\sum_{k=0}^{n-1}\{f(x_k)+f(x_{k+1})\}+\frac{b-a}{2}\sum_{k=0}^{n-1}f((x_k+x_{k+1})/2)\\=\frac{1}{2}T_n+\frac{b-a}{2}\sum_{k=0}^{n-1}f((x_k+x_{k+1})/2) T2n​=2b−a​k=0∑n−1​{2f(xk​)+f((xk​+xk+1​)/2)​+2f((xk​+xk+1​)/2)+f(xk+1​)​}=4b−a​k=0∑n−1​{f(xk​)+f(xk+1​)}+2b−a​k=0∑n−1​f((xk​+xk+1​)/2)=21​Tn​+2b−a​k=0∑n−1​f((xk​+xk+1​)/2)直观的说，新一次的迭代等于上次迭代值的一半加上新插值求和(再乘以$\frac{b-a}{2}$)。能统一公式的话，每次迭代就好办了，最后的判断条件：$$|T_n-T_{2n}|<ϵ$$满足之后，这有$${\int }_a^{b}f(x)dx=T_{2n}$$

变步长Simpson积分

辛普森积分可以由梯形积分来求解，推导省略【9】,根据上面的推导，假设辛普森面积积分用S表达，则有：$$S_n=\frac{4T_{2n}-T_n}{3}$$而$$S_{2n}=\frac{4T_{4n}-T_{2n}}{3}$$说明三次梯度结果能得出两次辛普森迭代结果，其他方法和梯度一致。

参考文档

1.GNU Scientific Library
2.Numerical Integration Using Trapezoidal Method C Program
3.C library for Numerical Integration
4.Numerical integration
5.Numerical Recipes:The Art of Scientific Computing
6.Simpson’s Rule （辛普森法则）
7.Integration
8.可视化微积分：一种全新的视野解释抛物线的面积
9.常用算法程序集(C语言描述)第三版