jiyang1971的个人博客分享 http://blog.sciencenet.cn/u/jiyang1971

博文

计算方法之常微分方程的初值问题 精选

已有 4125 次阅读 2018-12-3 11:36 |个人分类:大众物理学|系统分类:科普集锦

 


 

靡不有初,鲜克有终。

 

 

 

大学课本里的问题(不管是例题还是习题),大多是几十年乃至几百年前的问题。时代变了,解决问题的思想和原理可能还没变,但是相应的技术和方法却变了很多。比如说,力学只是求解二阶微分方程而已(牛顿第二运动定律),以前追求解析表达式、各种微扰展开,因为那时候不能做太多太繁琐的计算工作。随着计算机的发明和计算能力的突飞猛进,现在处理实际问题时主要依靠数值解法。所以,有必要对此有些了解,我们这个计算方法系列就再拿常微分方程的初值问题做个例子吧。

问题很简单,就是求解一阶微分方程$y'=\frac{dy}{dx}=f(x,y)$,实现知道函数$y$在某个位置处的初始数值$y(x_0)=y_0$。知道加速度(或者力)求速度,知道速度求位移,都可以归结为这种问题。

如果$f(x,y)$的形式很简单(比如,不显含$x$$y$,或者可以表达为$g(x)h(y)$的形式,或者通过变量代换变成这种简单的形式,等等),通常是可以求得解析表达式的。但是解析表达式的计算并不一定就简单,在很多情况下,甚至不可能求得$y$的解析解,就用一系列离散的数值$(x_i,y_i)$(其中$i=1,2,3......$)来近似地表示$y$。对于等间隔选取的$x_i$,就有$y_i=y(x_i)=y(x_0+ih)$

我们知道,在$x_i$处有泰勒公式,$y_{i+1}=y_i +y'(\xi)h$,其中$x_i<\xi<\x_{i+1}。关键在于找到$\xi$,这里有很多选择,也就有了很多不同的算法。

 

最简单粗暴的方法就是选择$\xi=x_i$,那么$y_{i+1}=y_i +y'(x_i)h=y_i +f(x_i,y_i)h$。注意到$x_0$处的初值和导数都是已知的(分别为$y_0$$f(x_0,y_0)$),自然就可以得到下一个位置的数值$y_1=y_0+f(x_0,y_0)h$,以及相应的导数$f(x_1,y_1)$(这就是由01)。接下来可以求出下一个位置,由12,由23,以至于无穷。这就是欧拉法。

也可以选择$\xi=x_{i+1}$,即$y_{i+1}=y_i +y'(x_{i+1})h=y_i +f(x_{i+1},y_{i+1})h$。这就稍微麻烦一些,因为导数依赖于待求的数值$y_{i+1}$,所以这时需要求解关于$y_{i+1}$的函数方程,但是可以做到的。这也是欧拉法,称为隐式的欧拉法,上面那个自然就是显式的欧拉法。

上面两种方法各有侧重,一个用了前端值,一个用了后端值,完全可以再公平一些嘛——采用二者的平均值来近似导数值,进而得到微分方程的解。这就是改进的欧拉法,具有预估校正系统的欧拉法。具体地说,先用前端值算出一个预估值$y_p=y_i+f(x_i,y_i)h$,再用这个预估值算出一个新值$y_c=y_i+f(x_{i+1},y_p)h$,最后取二者的平均值,作为下一个位置上的函数值$y_{i+1}=(y_p+y_c)/2$

还可以选择$\xi=(x_{i+1}+x_i)/2$,这相当于把分区加倍了。还是换种方式来说吧。考虑$y'(x_i)=(y(x_{i+1}-y(x_{i+1})/2h$,就可以得到$y_{i+1}=y_{i-1}+f(x_i,y_i)2h$。注意,这里需要附加条件,因为从0开始直接就蹦到2去了,必须先想办法得到1的值才行。以后的每次计算也都涉及此前两个位置的数值。这还是欧拉法,所谓的两步格式的欧拉法。这个方法当然也可以用预估矫正系统来改进。

 

如果函数$f(x,y)$具有更好的性质(欧拉法通常只要求函数连续,且满足所谓的Lipschitz条件),就可以更好地估计导数值,从而得到更好的算法——这就是龙格-库塔方法(Runge-Kutta)的精神。具体地说,令$y'(\xi)=\sum _i \lambda _i K_i$,加权系数$\lambda _i$满足归一化条件$\sum _i \lambda _i=1$,而$K_i$一些导数的预估值。有一套具体的方法来选择这些预估值和加权系数,使得计算结果的截断误差达到最小,任何标准的计算方法教材都可以找到,这里就不介绍了。这里只需要指出两点:

首先,经典的龙格-库塔方法是一种常用的四阶方法,具有很高的精度,其形式为

$y_{i+1}=y_i+(K_1+2K_2+2K_3+K_4)h/6$

$K_1=f(x_i,y_i)$

$K_2=f(x_i+h/2,y_i+K_1h/2)$

$K_3=f(x_i+h/2,y_i+K_2h/2)$

$K_4=f(x_i+h/2,y_i+K_3h)$

其次,龙格-库塔方法要求函数$f(x,y)$具有很好的光滑性,如果这个条件不能满足,其得到的结果甚至可能还不如改进的欧拉法。

 

龙格-库塔方法的精度很高,但是需要多次计算导数值,有些麻烦(就是说计算量大)。计算了$k$个数值后,如果能够把这些数值充分利用起来,第$k+1$个数值就可能估计得更精确——这就是线性多步法的精神。具体地说,就是用前面的导数值来估计当前的导数值,$y_{i+1}=y_i+[\sum _k \beta _k f(x_{i-k},y_{i-k})]h$。四阶显式(Adams格式)是$y_{i+1}=y_i+[55y'_i-59y'_{i-1}+37y'_{i-2}-9y'_{i-3}]h/24$

其精度为$h^5$。这种方法当然也有隐式和改进的格式(预估校正),还有其他变种,如Simpson格式、Milne格式和Hamming格式等,就不一一介绍了。

 

能解一阶微分方程,自然也就能解高阶微分方程,能解单变量的微分方程,自然也就能解多变量的微分方程(也就是偏微分方程)。只需要注意到,$n$阶微分方程可以变量代换表示为$n$个一阶微分方程(要点在于把$y'$视为与$y$无关的独立函数$z=y'$,更高阶的导数也是如此),而$n$变量的偏微分方程,只需要暂时把其他$n-1$个变量视为不变量,就成了单变量的微分方程。这样,问题就归结为$n$个一阶微分方程构成的微分方程组了。

至于说这些计算方法的适用条件和计算精度,都需要进一步的分析,难倒是不难,就是有些繁琐,随便哪本书里都有的,我就偷个懒吧。



http://blog.sciencenet.cn/blog-1319915-1149589.html

上一篇:一些诗词方面的书
下一篇:计算方法之引力波数据的初步显示(多图)

12 李颖业 康建 郑永军 刘安金 曾泳春 赵克勤 李毅伟 张江敏 徐磊 朱伯靖 林楠 苏德辰

该博文允许注册用户评论 请点击登录 评论 (11 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备14006957 )

GMT+8, 2018-12-19 23:39

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部