二阶偏微分方程组 龙格库塔法_数值积分公式的推广:龙格库塔法(一)

1. 写在前面的话

龙格库塔法是求解常微分方程的数值计算方法。它始于1895年,历史悠久,后人以龙格(Runge)和库塔(Kutta)两位计算数学家的名字为它命名。关于它的研究文献有很多,至今仍有新的论文不断涌现,卷帙之浩繁,恐怕难以完全统计。

先后在此领域内作出重要贡献的科学家除了龙格和库塔两位外,还有霍伊恩(Heun),奈斯特龙(Nyström),胡塔(Huta),吉尔(Gill),柯蒂斯(Curtis),布切尔(Butcher),默森(Merson),库珀(Cooper),海尔(Hairer),万纳(Wanner),桑兹塞尔纳(Sanz-Serna)等。

其中新西兰数学家布切尔是该问题的当代大师. 他创造性地将抽象代数和图论的方法引入到龙格库塔法的研究当中,奠定了龙格库塔法根树表示论的理论基础,进而改变了整个数值方法的面貌。以他的名字命名的数学名词有:布切尔群,布切尔表,布切尔级数(简称为B级数)。关于龙格库塔法,除了一系列重要的研究性论文,他还写过好几篇值得一读的综述类文章,并有专著传世,感兴趣的读者可以查阅其著作或文献(比如【1,2】)。

本文仅就龙格库塔法背后隐藏的机理作一番初步探讨,有兴趣的读者可以进一步思考或加以研究。

2. 数值积分公式

在介绍龙格库塔法之前,我们先来回顾一下数值积分公式。随便翻开一本计算方法的教科书,大家都能看到下面的数值求积公式:5a2bd5da97da47e5f264ac5399636095.png其中 941dc6063a1da534c53d0dcdbe2f4021.png 称为权系数,而 72a9361142d94d1308b2830b53cf9482.png 称为求积节点。这里假定被积函数比较光滑。通常情况下,节点在积分区间 (0, 1)内部或边界上选取。

其中最常用的求积公式是所谓的插值型求积公式,即用函数 f 的拉格朗日插值多项式来逼近被积函数, 从而得到a2eca231362995a612c0953596272a61.png再将积分号与求和符号交换, 得到5a2bd5da97da47e5f264ac5399636095.png其中873ea6111623c9e56afe759c9a6b2ca7.png注意, 当节点个数只有一个时 (即 s = 1),可定义 6af6d4768a919ee7eb771aa5638dd4c5.png

求积节点的选取可以多种多样,不同的选取方式将导致不同的逼近误差,最终影响到数值求积的精度。比如,若分别选取高斯(Gauss)节点、拉道(Radau)节点、洛巴托(Lobatto)节点,则相应得到高斯型、拉道型及洛巴托型数值求积公式,其代数精确度分别为 fdc9c4e5514a1a2ad4e23daa397211dc.png。所谓数值积分公式的代数精确度,是指当被积函数 为多项式函数时,使得数值求积公式精确成立的多项式的最高次数。

显然,该数值积分公式可以直接用于定积分的近似计算。对于简单的 维一阶常微分方程初值问题df810928e6b1faa3f883b5f3df793aa7.png

这里字母上方的点表示对时间变量 t求导,若需计算e6850bc732ef359fbb88b7fd3dee3274.png时刻的值,则可先通过积分运算得到解的解析表达式692c72be67f7edbfa56e5ad99b4f6cd6.png其中 097af4bb8d334f524ba9de61dc9cff2d.png. 然后直接使用数值积分公式求积分即可。

3. 数值积分带给我们的启示

现在让我们来考虑下面的一阶常微分方程初值问题60bede8027e37f72a8eed382a1aa1e23.png为了保证解的存在唯一性,通常假定右端函数足够光滑且关于z 满足李普希茨条件。

注意,此方程组不是简单的常微分方程组,右端函数 f 依赖于函数 z,但本身 未知,一般不能直接求出解析解 -- 当然不排除某些特殊方程组可以利用某些特别的技巧求得解析解,但大多数情况下是无法得到解析解的。

因此,实际当中我们往往会考虑使用合适的数值方法来求近似解。设精确解 z 在一系列离散时间点 ab20d9c783c4041926b0dc01cbc85cb4.png处的近似值为 6734bc235217f81c27b783e8cc1816f2.png,我们称由 a9d9774ec13621743ad042af25dec999.png 计算46fed86410cb28386367a12d1e71993b.png 的数值方法为单步法。下面我们考虑使用数值积分公式构造单步法。

首先将时间区间进行分割,得到离散时间区间 7a9802294eb091d590dc68b8331d5af0.png。 不失一般性,下面只考虑第一个时间区间的情况,设 097af4bb8d334f524ba9de61dc9cff2d.png为步长。通过对原微分方程两边求积分,代入初值条件,原方程可化为积分方程

e9f889a0b6ef080512fda613a11923b1.png

我们希望根据初值 求出下一步的值,为此, 在上式中令 026a12ce441199f5f0ae4457e30121dc.png,得到

c93391de83ac95db521d1d02750b1e9e.png

右边的积分很容易让我们联想到上一节的数值积分公式。特别地,如果我们使用左矩形数值积分公式(即只含有一个求积节点696923c071febfe3a6445c39b831baa8.png 的情况),则有

d35a877d47acb423fe12502d32a2196a.png

将精确解的符号替换成相应数值解的符号,并强行改成等式,即得显式欧拉法c12665b9438c7fbe5c84e71d49899593.png利用上式可由初值 50ab43f786de7bdec5358e431df8173c.png 求出 8095ba39348d3c7179ae0bd456bdfe6a.png,然后将 8095ba39348d3c7179ae0bd456bdfe6a.png作为下一步的初值,以此类推可以求得各个离散时间点处的近似值。显式欧拉法又称为欧拉折线法,是大数学家欧拉(Euler)于1768年提出来的,也是计算方法入门教科书上必提到的方法之一,通常还会有直观的图解(参见 ODE 数值方法(1):方法简介)。

遗憾的是,显式欧拉法只有一阶精度,稳定性也不太好。如果考虑使用右矩形数值积分公式,则可得到隐式欧拉法,稳定性比显式欧拉法要好很多,但精度也只有1阶。

科学研究总是从最简单的情况开始,然后试着考虑更为复杂的情况。一个很自然的想法是,为了提高精度,我们可否使用更多的求积节点?

一般地,我们在 [0, 1] 内选取 个求积节点,并利用相应的数值积分公式,得到ed51e97757b0ae3eaac48632bb422cd5.png

fc9eea336685c80dc42a6968229e390d.png依样画葫芦,改写后得5c6e81dd2d72801fb49ba0d65a2f4d5a.png其中7c8523bf90b0b5f181829332d957797b.png是积分点处函数值的近似。

问题来了,我们不知道 3758da22d931d429d4db615812937ae1.png 到底取何值,只是猜测它应该是真解 70494843063d779edb41c7391af61152.png的某种近似。幸运的是,我们有精确解的表达式。

对每个 8e7b5ec301187612225b0f1ca1801ac0.png8fdc1d6e58081ef7ad1fef062f07d9fb.png

接下来,我们需要想一想怎么把它跟 3758da22d931d429d4db615812937ae1.png联系起来,自然还是想到用数值积分公式。可是,上面的积分不是在 (0, 1) 上进行的,好像不能直接套用前面的数值积分公式。

有人可能会马上说,这还不简单,把积分区间通过变量替换变成  (0, 1) ,然后再使用数值积分公式不就行了吗?这里要注意一点,我们自然还是希望只用到原来的那些节点,而不希望引进新的额外的节点,否则问题复杂化。因此,稍微仔细一想,若按照上述替换积分变量的思路,使用数值积分后节点的位置必然会由于伸缩变换而发生错位,与原来的节点不一致,不但无法跟原来节点处的3758da22d931d429d4db615812937ae1.png 联系起来,还引入了新的节点处的未知量!

那么,接下来该怎么办呢?回想一下,插值型数值求积公式的本质是将被积函数替换成其拉格朗日插值多项式再积分!现在我们故技重施,把被积函数换成它基于积分点的插值:对 bb65c991a8e4908d7323e33654e77068.png

56bbc37317136dcca5eb2b535bd129da.png

615e75386cec8c4df5c773df85e42677.png节点错位的问题消除了!引进记号f89285be1b97460f382a8ebd2784ca04.png并将上式改写后得到739acf3791cf8da01c1ef0db41954dc5.png也就是说,我们虽然还是没有得到 3758da22d931d429d4db615812937ae1.png的显式表达式,但我们得到了关于 s个未知量 3758da22d931d429d4db615812937ae1.png s 个方程式,可以通过迭代法来数值地求出每个 3758da22d931d429d4db615812937ae1.png

综上所述, 我们得到下面的单步法:f65c893cf99e504f15496741604249b8.png

5c6e81dd2d72801fb49ba0d65a2f4d5a.png

其中 dab4e3d061e701a6a909a56e48add792.png 如前定义。为了得到下一步的数值解,我们需要先计算 s 个未知量 3758da22d931d429d4db615812937ae1.png,即求解前面含 s 个未知量的代数系统,然后代入最后一个等式求得 。这就是龙格库塔法。

可以证明,上述单步法的精度与相应的数值积分公式的精度是一致的。特别地, 如果我们选取求积节点为高斯(Gauss)节点 (见下图),则我们得到著名的高斯-勒让德(Legendre)法,它具有最高阶的局部截断误差阶 ad4822b7233121753f238073bafee92c.png

d7fb444c8be425325913f7cfa45ae302.png

前面构造的这个单步法,其实等价于所谓的配置法(也叫配点法),配置法的定义如下:设节点 053c472b40536bb649c16fad861f96c1.png为两两互不相同的实数。令 次多项式 f712c02d4f7630ea2c6f09296bcae000.png满足以下条件20f8b7dd7ec7643c795e7680019a00ab.png定义下一步的数值解为 f0652d96b80740912823289ad0be3bfa.png,则称此单步法为配置法。配置法的思想比较简单直观,它只要求在 s(有限)个点处满足原微分方程,而原微分方程是要求在整个时间区间(无穷多个点处)上精确成立。其实这里体现了计算数学中常见的一种逼近思想,用有穷逼近无穷,用有限代替无限。

请读者思考为什么上述两种方法是等价的。

(未完待续)

参考文献

【1】J.C. Butcher, G. Wanner, Runge-Kutta methods: some historical notes, Appl. Numer. Math.22 (1996), 113-151.


【2】J.C. Butcher, The Numerical Analysis of Ordinary Differential Equations: Runge–Kutta and General Linear Methods, John Wiley & Sons, 1987.


作者简介

7c2956970a77d1fb077c22b9fd159d6a.png

唐文生,长沙理工大学数学与统计学院, 讲师   

研究兴趣:哈密尔顿辛几何算法.

近年来,基于正交多项式展开技巧发展了连续级算法的较一般性构造理论,并将之用于动力系统保结构算法的研究,提出了构造几类几何积分子的全新方法. 首次发现利用数值通量的不同选取可以构造不同的保持辛结构的间断Galerkin方法,证明了方法的超收敛性,并建立了其与连续级算法的内在联系。主持并完成国家自然科学基金青年项目1项,湖南省教育厅项目1项。

  • 0
    点赞
  • 0
    评论
  • 1
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

打赏
文章很值,打赏犒劳作者一下
表情包
插入表情
评论将由博主筛选后显示,对所有人可见 | 还能输入1000个字符
相关推荐
©️2020 CSDN 皮肤主题: 数字20 设计师:CSDN官方博客 返回首页

打赏

Jack Weavi

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者