大学课本里的问题（不管是例题还是习题），大多是几十年乃至几百年前的问题。时代变了，解决问题的思想和原理可能还没变，但是相应的技术和方法却变了很多。比如说，力学只是求解二阶微分方程而已（牛顿第二运动定律），以前追求解析表达式、各种微扰展开，因为那时候不能做太多太繁琐的计算工作。随着计算机的发明和计算能力的突飞猛进，现在处理实际问题时主要依靠数值解法。所以，有必要对此有些了解，我们这个计算方法系列就再拿常微分方程的初值问题做个例子吧。
问题很简单，就是求解一阶微分方程y′=dydx=f(x,y)，实现知道函数y在某个位置处的初始数值y(x0)=y0。知道加速度（或者力）求速度，知道速度求位移，都可以归结为这种问题。
如果f(x,y)的形式很简单（比如，不显含x或y，或者可以表达为g(x)h(y)的形式，或者通过变量代换变成这种简单的形式，等等），通常是可以求得解析表达式的。但是解析表达式的计算并不一定就简单，在很多情况下，甚至不可能求得y的解析解，就用一系列离散的数值(xi,yi)（其中i=1,2,3......）来近似地表示y。对于等间隔选取的xi，就有yi=y(xi)=y(x0+ih)。
我们知道，在xi处有泰勒公式，yi+1=yi+y′(ξ)h，其中$x_i<><>
最简单粗暴的方法就是选择ξ=xi，那么yi+1=yi+y′(xi)h=yi+f(xi,yi)h。注意到x0处的初值和导数都是已知的（分别为y0和f(x0,y0)），自然就可以得到下一个位置的数值y1=y0+f(x0,y0)h，以及相应的导数f(x1,y1)（这就是由0而1）。接下来可以求出下一个位置，由1而2，由2而3，以至于无穷。这就是欧拉法。
也可以选择ξ=xi+1，即yi+1=yi+y′(xi+1)h=yi+f(xi+1,yi+1)h。这就稍微麻烦一些，因为导数依赖于待求的数值yi+1，所以这时需要求解关于yi+1的函数方程，但是可以做到的。这也是欧拉法，称为隐式的欧拉法，上面那个自然就是显式的欧拉法。
上面两种方法各有侧重，一个用了前端值，一个用了后端值，完全可以再公平一些嘛——采用二者的平均值来近似导数值，进而得到微分方程的解。这就是改进的欧拉法，具有预估校正系统的欧拉法。具体地说，先用前端值算出一个预估值yp=yi+f(xi,yi)h，再用这个预估值算出一个新值yc=yi+f(xi+1,yp)h，最后取二者的平均值，作为下一个位置上的函数值yi+1=(yp+yc)/2。
还可以选择ξ=(xi+1+xi)/2，这相当于把分区加倍了。还是换种方式来说吧。考虑y′(xi)=(y(xi+1−y(xi+1)/2h，就可以得到yi+1=yi−1+f(xi,yi)2h。注意，这里需要附加条件，因为从0开始直接就蹦到2去了，必须先想办法得到1的值才行。以后的每次计算也都涉及此前两个位置的数值。这还是欧拉法，所谓的两步格式的欧拉法。这个方法当然也可以用预估矫正系统来改进。
如果函数f(x,y)具有更好的性质（欧拉法通常只要求函数连续，且满足所谓的Lipschitz条件），就可以更好地估计导数值，从而得到更好的算法——这就是龙格-库塔方法（Runge-Kutta）的精神。具体地说，令y′(ξ)=∑iλiKi，加权系数λi满足归一化条件∑iλi=1，而Ki一些导数的预估值。有一套具体的方法来选择这些预估值和加权系数，使得计算结果的截断误差达到最小，任何标准的计算方法教材都可以找到，这里就不介绍了。这里只需要指出两点：
首先，经典的龙格-库塔方法是一种常用的四阶方法，具有很高的精度，其形式为yi+1=yi+(K1+2K2+2K3+K4)h/6
K1=f(xi,yi)
K2=f(xi+h/2,yi+K1h/2)
K3=f(xi+h/2,yi+K2h/2)
K4=f(xi+h/2,yi+K3h)
其次，龙格-库塔方法要求函数f(x,y)具有很好的光滑性，如果这个条件不能满足，其得到的结果甚至可能还不如改进的欧拉法。
龙格-库塔方法的精度很高，但是需要多次计算导数值，有些麻烦（就是说计算量大）。计算了k个数值后，如果能够把这些数值充分利用起来，第k+1个数值就可能估计得更精确——这就是线性多步法的精神。具体地说，就是用前面的导数值来估计当前的导数值，yi+1=yi+[∑kβkf(xi−k,yi−k)]h。四阶显式（Adams格式）是yi+1=yi+[55y′i−59y′i−1+37y′i−2−9y′i−3]h/24，其精度为h5。这种方法当然也有隐式和改进的格式（预估校正），还有其他变种，如Simpson格式、Milne格式和Hamming格式等，就不一一介绍了。
能解一阶微分方程，自然也就能解高阶微分方程，能解单变量的微分方程，自然也就能解多变量的微分方程（也就是偏微分方程）。只需要注意到，n阶微分方程可以变量代换表示为n个一阶微分方程（要点在于把y′视为与y无关的独立函数z=y′，更高阶的导数也是如此），而n变量的偏微分方程，只需要暂时把其他n−1个变量视为不变量，就成了单变量的微分方程。这样，问题就归结为n
个一阶微分方程构成的微分方程组了。
至于说这些计算方法的适用条件和计算精度，都需要进一步的分析，难倒是不难，就是有些繁琐，随便哪本书里都有的，我就偷个懒吧。
评论 (0)