数值积分，暴力求解微分方程
在学习高数的过程中，相信很多人都“清楚”地认识到大部分积分式不可求解 （准确来说是懒得解），积分式内的结构就是一个dy/dx = f(x)的最基本的微分方程。
但是现实是大部分一阶微分方程都是dy/dx = f(x, y)的结果，也就是说目标方程本身就在其导数内，such：dy/dx=y，自觉上(查表)知道y=αe^t，但是如果真的要求解的话也会让很多人无从下手。虽然有各种求解微分方程的办法，但是都异常繁琐且不通用。这时候就到了数值积分发挥用处的时候了。
什么是数值积分
为了更贴合使用数值积分的情况，我们来想象一个情景：一辆汽车在直路上行走，路上有一个原点，并且汽车的速度与汽车到原点的距离成正比，在开始时(t=0)，汽车距离原点1米，求汽车与原点的距离y与时间t的关系。
根据题意不难列出：v=dy/dt=αy, y|t=0 = 1，当然也不难知道y(t)=e^αt，但是假如我们特别蠢，完全不会求解这个方程，但是又特别想知道结果是怎么样的，那么就只能暴力计算：算出当前速度，把汽车按照这个速度前进一点点距离，然后再计算速度，然后再前进一点点，然后......
这种做法就是数值积分的核心思想了，数值积分的核心思路就是积分的原理:
这时候就会有人说，每一步都会累积误差，完全不是精确的，怎么能叫求解呢。问得好，尽管是有误差，但是精度可以控制在接受范围内，科学家们已经用这种办法对几百亿个星体进行演算模拟出星系碰撞甚至宇宙的演化。
了解了数值积分是什么后，就可以开始了解求解细节和优化的办法了
基础 -- 欧拉法
虽然第一个是欧拉发明的方法，但却是最简单的一个办法
注意：以下所有方程都为 dy/dt=f(t,y) 的结构，且附带初值条件y0=○。目标方程y(t)可以为向量函数，也就是说所有方法都适用于高维空间
看到导数的定义式：
可以得到：
我们记y_n为当前知道的点，对应的时间为t_n，y_n+1为下一个要求的点，那么算法可以写为：
进阶 -- 中点法
有些人可能看出来了一个问题：如果遇到曲率很大的地方，直接采用y_n点的导数会与实际曲线存在很大偏差
对此，我们又想到了一样东西：一段圆弧的中点的切向量总与两端点连线的向量平行
于是得到中点法方法:
对于欧拉法和中点法分别对上面汽车的例子作图来看看
可以中点法比欧拉法更接近目标方程, 但是也有比较大的误差.
接下来的方法已经没有明显的对比效果了，在最后我会用“双蝴蝶”做一个总结
究极 -- Heun法（梯形法）
看到不知道什么什么定理（我是真的忘了名字），F是f积分后的原函数
因为不能确定ξ的值，所以一般用ab两点f的平均数代替f(ξ)，即：
根据这个思路得到：
这个方法对一种叫做刚性系统的东西有特攻，后面再详细讲述, 并且善于利用的话可以做到任意小的误差
我也不知道是什么原理的方法 -- Runge-Kutta法（rk4）
rk4的计算方法如下：k1是y_n起点的导数，k2的第一中点的导数，k3是第二中点的导数，k4是k3终点的导数，各自加权后得到最终斜率
欧拉法的误差为Δt^2，梯形法的误差为Δt^3，而rk4的误差达到了Δt^5
因为计算量少，并且误差小，所以rk4是一栋非常常用的求解微分方程的方法
Butcher tableau 和 各种显式rk方法
受rk4的启发，人们开始拓展多点权重累加的方法，通式如下：
为了简洁描述，发明了一种叫Butcher tableau的表格：
有了这种表格可以自己尝试调整系数, 去寻找快捷方便的求解方法. 但是为了保证方法是自洽的(也就是不会自身出问题), 需要有以下条件:
简单来说就是：每一行的a相加等于c，且b相加等于1
以下为各种不同的显式rk方法 en.m.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods
上述方法在大部分情况都表现不错，但是缺点是需要手动调整时间间隔以确保性能精度兼得。以懒著称的数学家发明了根据误差自动调整间隔的自适应方法
误差e的计算就是用不同的权重计算y_n+1，如果两个结果相差太大则减少间隔以确保精度。
这里我需要说一声抱歉，我一直都找不到如何根据误差调整间隔，如果有读者知道的话请务必在评论区大概讲一下是怎么工作的
什么是刚性系统
在某些条件下，使用自适应方法求解时，会出现间隔不断减少，以致需要花费巨量算力才能得出结果，或者根本得不出结果。这种系统就叫做刚性系统。刚性与目标方程的形状无关，而与系统的特征值有关，特征值的实部越靠近负无穷刚性就越大。通常来说特征值的实部小于-2就为刚性系统。如何计算系统特征值？问得好，我也不知道
🌰：dy/dt=-15y就是一个典型的刚性系统，当t=0时，y=1
分别使用间隔为0.25和0.125的欧拉法求得的图像：
可以看到没有收敛在目标函数附近，而使用间隔为0.125的rk4方法可以得到:
如此可以看到部分方法对刚性系统还是有效果的
针对刚性系统，rk党并没有放弃治疗，并且研发出了
隐式rk方法
显式rk方法是从上往下计算，非常简洁方便。而隐式rk方法是每条式子都包含自身和其他式子
尽管隐式rk法可以普遍地求解刚性和非刚性系统，但是相比显式rk方法，隐式需要的算法复杂度和算力真的大太多了。需要了解如何计算隐式rk方法的可以去查阅相关文献
各种隐式rk方法：
预测-评估-校正-评估方法
实际上梯形法就是最简单的预测评估-校正-评估方法
步骤如下: 1.预测: 使用欧拉法预测下一个点 2.评估: 计算下一个点的斜率 3.校正: 欧拉法的斜率与评估的斜率作平均得到新斜率 4. 评估: 计算新斜率对应下一个点的斜率 (一般来说不需要最后这一步)
这就是PECE的步骤了, 当然最后评估得到的斜率也可以继续拿回去纠正, 这样就变成了PECECE方法
除非在刚性非常大的系统内, 这种方法会随着纠正的次数增多而逐渐接近目标方程. 当做了无穷部评估校正步骤后, 必定会落在目标方程上
应用上述方法也是非常简单的, 如图:
分别使用欧拉法, 梯形法, 和rk4方法求解著名的洛伦兹方程: (间隔: .001, 起点:[.1, .1, .1])
红色 - 欧拉法 绿色 - 梯形法 蓝色 - rk4方法
下一个项目应该就是随机抽取一种幸运方法拿去模拟三体运动了
彩蛋：一种简单粗暴地把数据平滑化的方法
无论在做什么项目，比如说上述的数值积分，总会得到一堆离散的数据，有时候就是不想用线性插值，想要各个数据直接平滑过渡，那么这里有一种懒人方法：多项式
如果现在有n组数据（这里用[ti, yi]做例子），那么假设光滑的目标函数为一元n-1次方程：
那么可以列一个有n个方程度方程组，求出各系数就得到了完全光滑的目标函数y(t)，不过需要注意，这个函数只有在最小t和最大t之间才好看，超过了这个范围将会急速变为无穷。
什么？这么大的方程组不会解？真的没办法呢，列一个这样的线性方程组：
然后让计算机解这个就行啦
