欧拉方法

维基百科,自由的百科全书
跳转至: 导航搜索

数学计算机科学中,欧拉方法,命名自它的发明者萊昂哈德·歐拉,是一种一阶数值方法,用以对给定初值的常微分方程(即初值問題)求解。它是一种解决数值常微分方程的最基本的一类显型方法(Explicit method)。

非正式的几何描述[编辑]

欧拉方法的图示。待求的曲线为蓝色,它的折線近似为红色。

考虑计算這樣的一个未知曲線的形状:它具有给定的起点并且满足一个给定的微分方程。 这里,所谓“微分方程”可以看作能够通过曲线上任意点的位置而计算出这一点的切线斜率的公式。

思路是,一开始只知道曲線的起点(假设为A_0),曲線其他部份是未知的,不過通过微分方程,A_0的斜率可以被计算出来,也就得到了切线。

顺着切线向前走一小步到点A_1。如果我们假设A_1是曲线上的一点(实际上通常不是),那么同样的道理就可以确定下一条切线,依此类推。在经过几步之后,一条折线A_0A_1A_2A_3\dots就被计算出来了。一般情况下,这条折线与原先的未知曲线偏离不远,并且任意小的误差都可以通过减少步长来得到(虽然对于刚性方程而言会比较复杂)。

欧拉方法的推导[编辑]

图示为方程y'=y, y(0)=1的数值积分。蓝色为欧拉法,绿色为中点法,红色为精确解y=e^t。所用步长為h=1.0
图示为同一个方程在步长h=0.25时的结果。可以看出中点法比欧拉法收敛更快。

以以下微分方程為例

y'(t) = f(t,y(t)), \qquad y(t_0)=y_0,

希望用 y 在點 (t0,y(t0)) 附近的線性近似來得到其近似解(也就是 y泰勒展開式的前二項)。利用時間 tn 時的數值,若用單步的欧拉方法,可得到時間 tn+1 = tn + h 時的近似值如下:

 y_{n+1} = y_n + hf(t_n,y_n).  \qquad \qquad

欧拉方法是一種顯型方法,也就是說 y_{n+1} 的解是 y_i, i \leq n 的顯函數。

欧拉方法可以求解一階的微分方程,而任何N階的微分方程都可以表示成一階的微分方程。

對於微分方程

 y^{(N)}(t) = f(t, y(t), y'(t), \ldots, y^{(N-1)}(t))

可以通過新設輔助變量 z_1(t)=y(t), z_2(t)=y'(t),\ldots, z_N(t)=y^{(N-1)}(t),得到以下的等價方程

 \mathbf{z}'(t)
  = \begin{pmatrix} z_1'(t)\\ \vdots\\ z_{N-1}'(t)\\ z_N'(t) \end{pmatrix}
  = \begin{pmatrix} y'(t)\\ \vdots\\ y^{(N-1)}(t)\\ y^{(N)}(t) \end{pmatrix}
  = \begin{pmatrix} z_2(t)\\ \vdots\\ z_N(t)\\ f(t,z_1(t),\ldots,z_N(t)) \end{pmatrix}

這是一個以\mathbf{z}(t)為變量的一階系統,因此可以用歐拉法求解,也可以使用其他的一階數值方法。[1]

应用例题[编辑]

设微分方程为 y'=y ,初始值为 y(0)=1,试用欧拉方法求 y_3的近似值,步长为 h=1

欧拉法為:

 y_{n+1} = y_n + hf(t_n,y_n).

首先求f(t_0, y_0)(当n = 0),f的定義為f(t,y)=y,因此有

 f(t_0,y_0)=f(0,1)=1.

透過以上步驟,求得解曲線在点(0,1)的切线斜率。回顾直線斜率的定义:y变化量和t变化量的比值,亦記作 \Delta y/ \Delta t

接著是

 h \cdot f(y_0) = 1 \cdot 1 = 1.
 y_0 + hf(y_0) = y_1 = 1 + 1 \cdot 1 = 2.

重复以上步骤求出 y_2y_3 的值。

 y_2 = y_1 + hf(y_1) = 2 + 1 \cdot 2 = 4
 y_3 = y_2 + hf(y_2) = 4 + 1 \cdot 4 = 8

由于欧拉法属于递归算法,把運算整理成表格也許有助於避免計算錯誤。

y_n t_n y'(t) h dy y_{n+1}
1 0 1 1 1 2
2 1 2 1 2 4
4 2 4 1 4 8

局部截尾误差[编辑]

欧拉法的局部截尾误差(Local truncation error, LTE)是指在实施一次欧拉法所产生的误差,是指经过一步的数值解y_1与在t_1 = t_0+h时精确解的误差。数值解y_1由以下给出:

 y_1 = y_0 + h f(t_0, y_0). \quad

对于精确解,我们使用泰勒级数展开给出:

 y(t_0 + h) = y(t_0) + h y'(t_0) + \frac{1}{2}h^2 y''(t_0) + O(h^3).

欧拉法的局部截尾误差为:

 \mathrm{LTE} = y(t_0 + h) - y_1 = \frac{1}{2} h^2 y''(t_0) + O(h^3).

y拥有三阶有界导数时,这个结果是成立的。[2]

结果显示:当步长h很小时,局部截尾误差近似与 h^2 成比例。也就是说,欧拉法没有其他的高阶方法如龙格-库塔法线性多步法精确,这些方法的局部截尾误差与h^pp>2)成比例。

全局截尾误差[编辑]

全局截尾误差(Global truncation error, GTE)是指在一个固定时间t时的误差,但是很多步之后该方法需要以从初始时间到达该时间来计算。全局截尾误差可以看做是一个每一步的局部截尾误差的累积效应。[3] 经过的步骤數為(t-t_0)/h,而每步的误差则正比于h^2。因此,可以预期全局截尾误差是正比于h的。[4]

这个直观的推测可以被嚴謹地證明。如果解y存在二阶有界导数,并且f關於y利普希茨连续的,那么全局截尾误差是有界的:

 |\text{GTE}| \le \frac{hM}{2L}(e^{L(t-t_0)}-1) \qquad \qquad

其中 M 是在给定区间内y的二阶导数的上界,Lf的利普希茨常数。[5]

这种精确的形式其实是没有什么意义的,通常情况下这个上界都會嚴重高估了欧拉法所造成的实际误差。[6]重要的是,这顯示了全局截尾误差是近似正比于h的,所以欧拉法被稱为是一阶的。[7]

註腳[编辑]

  1. ^ Butcher 2003,第3页; Hairer,Nørsett & Wanner(1993),第2页
  2. ^ Butcher 2003,第60页
  3. ^ Atkinson 1989,第344页
  4. ^ Butcher 2003,第49页
  5. ^ Atkinson 1989,第346页; Lakoba 2012,公式 (1.16)
  6. ^ Iserles 1996,第7页
  7. ^ Butcher 2003,第63页

参考文献[编辑]