[数值分析]欧拉格式和龙格库塔算法

发布于 2022-10-31  635 次阅读


首先明确欧拉格式和龙格库塔算法是用来干啥的,他们都用来解决同一个问题:

求解常微分方程

本文所作的的各种概念均有简化,可能与学术界不同,尽可能简单的理解。

常微分方程一般长这样:

图1.常微分方程的形式

可见y的导数由x和y共同表示。

现在我们的任务是求出y函数,或者说能够通过y0得到许多个y函数上的点y1, y2, y3 ... yn,从而可以用插值等方法构造出y的多项式来,总是就是要能够得到y的信息。

朴素欧拉格式(显式欧拉)

显式欧拉十分简单明了,先确定一个步长h,既然知道了x0,y0以及该点的导数值,那就假设这个y是一个完全线性的函数嘛,就像在直线上做的一样,不难得出以下公式:

图2.显式欧拉格式

只需要把k代入回去,就和书本上写的一样了。

但是这样得到的y[n + 1]显然十分粗糙,我们处理的大部分的函数并非线性。

隐式欧拉格式

我们将图2的公式公式整合一下得到以下的显式欧拉格式:

图3.显式欧拉格式

这里用的是x[n]和y[n]来得到y[n + 1],我们用的斜率是x[n]点的斜率,但是是不是也可以用x[n + 1]点的斜率呢?来试试吧。

改用x[n + 1]点的斜率写的欧拉格式如下:

图4.隐式欧拉格式

这是隐式欧拉格式,因为没有办法直接从右边算出y[n + 1],因为右边也有y[n + 1]。

怎么办呢?我们总得计算吧。

两步欧拉格式(优化一下!!)

那么就先用显示欧拉格式算出y[n + 1],再用这一个y[n + 1]去计算新的y[n + 1]。

写成公式是这样:

图5. 预估-校正系统的欧拉格式

这个就是预估-校正系统,其实也很好理解,先用显示欧拉预估一个y[n + 1],再用它去代入到隐式欧拉格式里计算y[n + 1],自己算自己属于是。

改进欧拉格式

前面的显式欧拉,隐式欧拉终究是线性的,代数精度最大也就1阶。

为了提高精度,我们需要对欧拉格式进行改进,回到常微分方程这里,我们如果对式子两边同时求x[n]到x[n + 1]的积分,可以得到:

图6. 改进欧拉格式

这就是利用梯形公式改进的欧拉格式。它具有2阶代数精度。

改进欧拉格式的截断误差和代数精度

首先看一下欧拉格式的截断误差定义:

图7. 截断误差定义

我们把y[n]的值代入进去,得到下面这个式子(为了方便计算,将n都加一):

我们将y[n + 1]和hf(x[n + 1], y[n + 1])分别进行泰勒展开得到如下式子,注意h * O(h^2) = O(h^3):

图8. 截断误差计算过程

所以改进欧拉格式具有2阶代数精度。

龙格库塔算法

龙格库塔不满足于两阶精度,他们想从精度的定义入手,得到更高精度的公式。

没错,从精度的定义入手十分关键。

类比高斯求积公式,我们前面用于改进欧拉格式的梯形公式也可以变形,变为对多个点进行加权平均的形式。

我们前面证明改进欧拉的泰勒展开最多展开到了O(h^3),如果我们展开到O(h^4)并且选出合适的参数使得误差值成立呢?

这就是龙格库塔的想法。

同时还有另外一个思路,我们要从y[n]求到y[n + 1],其实就是要确定一个比较好的斜率,使得求出的y[n + 1]尽可能准确,那么这个比较好的斜率就需要在x[n]到x[n + 1]中尽可能多的找点来做加权平均,将几个斜率加权平均得到较好的斜率,再从y[n]直接转移到y[n + 1]。

二阶龙格库塔(有很多种,这里写出其中一种,也就是改进欧拉格式):

三阶龙格库塔:

四阶龙格库塔:

感谢阅读。