跳去內容

歐拉法

出自維基百科,自由嘅百科全書
歐拉法嘅圖示。未知曲線為藍色,多邊形近似為紅色。

數學計算科學入便歐拉法(亦都叫做正向歐拉法)係一種一階數值過程,攞嚟求解啲具有畀定初始值嘅常微分方程(ODE)嘅。歐拉法係對常微分方程進行數值積分嘅最基本嘅顯式方法,亦都係最簡單嘅Runge-Kutta法。歐拉法以歐拉個名命名,佢喺佢本《計算積分學》(1768-1870年出版)一書入便有處理到。 [1]

歐拉法係一階方法,代表佢個局部誤差(每𨂾誤差)戥𨂾長嘅平方成正比、全局誤差(畀定時間嘅誤差)戥𨂾長成正比。歐拉法通常作為構造啲更複雜方法嘅基礎,譬如預測器-校正器法

非正式嘅幾何描述

[編輯]

考慮計算某條未知曲線個形狀;佢從畀定點開始、滿足一條畀定微分方程。基於呢層,微分方程即係一種公式,一旦計出任意一粒點個位置,透過呢條公式就計得出條曲線上嘅切線斜率嘅。

個惗法係,雖然條曲線初初都係未知嘅,但佢個起點係已知嘅(請參見右上方幅圖)。然之後,從微分方程就計得出條曲線喺度嘅斜率同埋切線。

沿住切線向上延伸多一點。沿住呢一細𨂾,個斜率唔會變化太大,所以好接近條曲線。若果假設個仲喺曲線上,戥上述喺點相同嘅推理過程都繼續使得。幾𨂾之後,就計出多邊形曲線。好多時,呢條曲線戥原始未知曲線嘅走差唔會太大;若果𨂾長夠細而且計算間隔有限,噉可以令到兩條曲線之間嘅誤差細啲: [2]

揀一個值對於每一𨂾同埋設置 。噉時,從係: [3]

ODE解嘅近似值,喺 嗰陣。 Euler法係顯式嘅,即對於,解嘅顯式函數。

就算Euler法集成唨一階ODE,但係任何N階嘅ODE都可以表示為啲一階ODE嘅系統。

爲唨處理以下方程:

可以引入輔助變量並得到等式:

呢個係變量入便嘅一階系統,可以用Euler法處理,或者實質上可以用任一其他啲都用於一階系統嘅方案處理。 [4]

例子

[編輯]

畀定初始值問題

諗住用歐拉法來逼近[5]

令𨂾長等於1(h = 1)

[編輯]
方程嘅數值積分圖示藍色係歐拉法;綠色中點法;紅色確切解, 𨂾長為h = 1.0 。

歐拉法係

所以首先必須計算 。喺條簡單嘅微分方程入便,函數由定義 。有初始值:

透過執行以上步驟,就搵到條切線嘅斜率、喺點戥方程解條曲線相切嘅 。考慮返斜率定義為變化量除以變化量 ,或者

下一步係將上述值乘上𨂾長 ,呢度等於1:

由於𨂾長係變化量,將𨂾長戥切線嘅斜率相乘嗰陣,就得到個變化值。然之後,呢個值會加到初始值以獲得下一個要攞嚟計算嘅值。

重複上述步驟即可搵到

由於呢個算法具有重複性,所以攞圖表形式嚟組織返啲計算會好有幫助,如下所示,嚟避免出錯。

0 1 0 1 1 1 2
1 2 1 2 2 4
2 4 2 4 4 8
3 8 3 8 8 16

個計算嘅結論係 。微分方程嘅精確解係 , 所以 。就算喺呢種𨂾長值大嘅特定情況之下,歐拉法嘅逼近唔點精確,之不過如圖所示,演算法個表現喺定性上係啱嘅。

MATLAB代碼示例

[編輯]
clear; clc; close('all');
y0 = 1;
t0 = 0;
h = 1; % try: h = 0.01
tn = 4; % equal to: t0 + h*n, with n the number of steps
[t, y] = Euler(t0, y0, h, tn);
plot(t, y, 'b');

% exact solution (y = e^t):
tt = (t0:0.001:tn);
yy = exp(tt);
hold('on');
plot(tt, yy, 'r');
hold('off');
legend('Euler', 'Exact');

function [t, y] = Euler(t0, y0, h, tn)
  fprintf('%10s%10s%10s%15s\n', 'i', 'yi', 'ti', 'f(yi,ti)');
  fprintf('%10d%+10.2f%+10.2f%+15.2f\n', 0, y0, t0, f(y0, t0));
  t = (t0:h:tn)';
  y = zeros(size(t));
  y(1) = y0;
  for i = 1:1:length(t) - 1
    y(i + 1) = y(i) + h * f(y(i), t(i));
    fprintf('%10d%+10.2f%+10.2f%+15.2f\n', i, y(i + 1), t(i + 1), f(y(i + 1), t(i + 1)));
  end
end

% in this case, f(y,t) = f(y)
function dydt = f(y, t)
  dydt = y;
end

% OUTPUT:
%     i    yi    ti    f(yi,ti)
%     0   +1.00   +0.00     +1.00
%     1   +2.00   +1.00     +2.00
%     2   +4.00   +2.00     +4.00
%     3   +8.00   +3.00     +8.00
%     4  +16.00   +4.00     +16.00
% NOTE: Code also outputs a comparison plot

R代碼示例

[編輯]

下面係使R編程語言編寫嘅示例代碼。

# ============
# SOLUTION to
#   y' = y,   where y' = f(t,y)
# then:
f <- function(ti,y) y

# INITIAL VALUES:
t0 <- 0
y0 <- 1
h  <- 1
tn <- 4

# Euler's method: function definition
Euler <- function(t0, y0, h, tn, dy.dt) {
  # dy.dt: derivative function
  
  # t sequence:
  tt <- seq(t0, tn, by=h)
  # table with as many rows as tt elements:
  tbl <- data.frame(ti=tt)
  tbl$yi <- y0 # Initializes yi with y0
  tbl$Dy.dt[1] <- dy.dt(tbl$ti[1],y0) # derivative
  for (i in 2:nrow(tbl)) {
    tbl$yi[i] <- tbl$yi[i-1] + h*tbl$Dy.dt[i-1]
    # For next iteration:
    tbl$Dy.dt[i] <- dy.dt(tbl$ti[i],tbl$yi[i])
  }
  return(tbl)
}

# Euler's method: function application
r <- Euler(t0, y0, h, tn, f)
rownames(r) <- 0:(nrow(r)-1) # to coincide with index n

# Exact solution for this case: y = exp(t)
#       added as an additional column to r
r$y <- exp(r$ti)

# TABLE with results:
print(r)

plot(r$ti, r$y, type="l", col="red", lwd=2)
lines(r$ti, r$yi, col="blue", lwd=2)
grid(col="black")
legend("top",  legend = c("Exact", "Euler"), lwd=2,  col = c("red", "blue"))

# OUTPUT:
#
#   ti yi Dy.dt         y
# 0  0  1     1  1.000000
# 1  1  2     2  2.718282
# 2  2  4     4  7.389056
# 3  3  8     8 20.085537
# 4  4 16    16 54.598150
# NOTE: Code also outputs a comparison plot

使埋其他𨂾長

[編輯]
相同插圖,h=0.25。

好似引言入便所述,若果𨂾長細啲,噉歐拉法亦都準確啲。下表顯示唨唔同𨂾長嘅結果。第一行對應於上一節入便個示例,第二行如圖所示。

𨂾長 歐拉法結果 誤差
1 16.00 38.60
0.25 35.53 19.07
0.1 45.26 9.34
0.05 49.56 5.04
0.025 51.98 2.62
0.0125 53.26 1.34

表格最後一欄記錄嘅誤差係精確解之間嘅差佮歐拉近似。喺份表格㞘底,𨂾長係上一行𨂾長嘅一半,而且啲誤差亦都係前一行啲誤差嘅一半咁上下, 表明誤差大致戥𨂾長成正比,至少對於相當細嘅𨂾長值嚟講。呢個對於其他方程式一般亦都係啱嘅。更多詳細信息請參睇全局截尾錯誤嗰部分。

其他啲方法、譬如圖入便亦都顯示有嘅中點法,表現得仲好:中點法嘅整體誤差戥𨂾長嘅平方大致成比例所以,歐拉法叫做一階方法,而中點法叫做二階方法。

由上表推得,要獲得啱嘅小數點後三位答案,所需嘅𨂾長約為0.00001,代表需要到400,000𨂾數。咁大量嘅𨂾數需要好高嘅計算成本。所以,啲人好多時採用高階嘅方法作爲替代,譬如啲Runge–Kutta法抑係啲線性多步法,尤其係要到高精度嗰陣。 [6]

推導

[編輯]

歐拉法可以透過多種方式推出。首先有上便啲幾何描述。

另一種可能性係考慮函數週邊嘅泰勒展開

條微分方程表明 。若果捉條式替換成泰勒展開式,而且忽略埋二次項同埋啲高階項,就得到Euler法。 [7]泰勒展開式可以似下低噉着攞嚟分析啲歐拉法產生嘅誤差,而且可以擴展到生成Runge–Kutta法

密切相關嘅推導係攞前向有限差分公式代替導數,

喺微分方程入便 。同樣,噉樣都推得出歐拉法。 [8]相似嘅計算推得出中點法同埋後褪歐拉法

最後,可以積分條微分方程式由,並應用微積分基本定理得到:

再透過左手矩形法(只用到一個矩形)來近似積分:

結合埋嗰兩個方程,又得到歐拉法。 [9]可以繼續採用呢種思路嚟得出各種線性多步法

局部截尾誤差

[編輯]

歐拉法嘅局部截尾誤差(LTE)係行親一𨂾造成嘅誤差,即係一𨂾之後數值解嗰陣嘅確切解之間嘅差。數值解由下式畀出

對於確切嘅解,使返上便嘅#推導部分入便提到嘅泰勒展開式:

歐拉法引入嘅局部截尾誤差由以下等式之間嘅差異畀出:

呢個結果有效,若果有有界三階導數。[10]

噉樣表明對於細 ,局部截尾誤差正比於 ,即會令到Euler法對於細啲嘅準確性降低,而唔似Runge-Kutta法同埋線性多步法之類其他啲高階竅,因為後者啲局部截尾誤差戥𨂾長嘅高啲嘅冪成正比。

一條略有唔同嘅局部截尾誤差公式,可以透過對於泰勒定理入便啲其餘項、用返拉格朗日形式獲得。若果有一個連續嘅二階導數,噉就存在一個令到:

[11]

喺上便嘅誤差表達式裏便,未知精確解嘅二階導數可以用一條表達式、涉及微分方程右側嘅嚟代替。根據方程有:

[12]

全局截尾誤差

[編輯]

全局截尾誤差(GTE)係固定時間嘅誤差 ,但係係經過好多𨂾步後,先從初始時間嗰度開始計。全局截尾誤差係每一𨂾入便整出嘅嗰啲局部截尾誤差嘅累積影響。 [13]𨂾數好容易確定為 ,之正比於 ,而且每一𨂾裏便整出嘅誤差正比於 (請參閱上一節)。所以,預期當中,全局截尾誤差會正比於[14]

呢種直觀嘅推理可以精確化。若果個解有有界二階導數、若果喺佢入便個第二個參數係Lipschitz連續嘅,噉個全局截尾誤差由條式得出:

其中係喺畀定間隔下、個二階導數嘅上限,而嘅Lipschitz常數 。 [15]

呢個界限嘅精確形式幾乎冇乜實際意義,因為喺大多數情況之下,界限遠遠高估唨Euler法造成嘅實際誤差。 [16]重要嘅衹係佢表明到個全局截尾誤差(大約)正比於 。所以,歐拉法着認為係一階嘅。 [17]

數值穩定性

[編輯]
用歐拉法計算嘅𨂾長 (藍色方塊)同埋 (紅色圓圈)。黑色曲線顯示唨確切嘅解。

Euler法喺數值上亦都可能係唔穩定嘅,尤其係對於剛性方程嚟講;噉樣即係對於冇精確求解嘅方程,數值解會變得非常之大。呢層可以用返線性方程來說明:

確切解係 ,隨住 而到零。但係,若果歐拉法以𨂾長應用於呢條式,噉個數值解喺定性上都係錯嘅:個數值解會振盪並增長(見圖),即唔穩定噉解。又如,若果𨂾長細啲、 ,噉個數值解又確實會衰減為零。

粉紅色嘅圓盤顯示唨Euler法嘅穩定區域。

若果將歐拉法應用於線性方程 ,噉呢個乘積嘅數值解會唔穩定、若果乘積喺以下區域之外:

似右圖所示,嗰柵區域叫做(線性)穩定性區域[18]喺呢個例子入便係-2.3,所以若果然之後,喺穩定區域之外,所以數值解唔穩定。

考慮到呢種限制與及𨂾長h誤差收斂之慢,代表除唨作為數值積分嘅簡單示例之外,Euler法唔經常會使到。

捨入誤差

[編輯]

到目前為止嘅討論都忽略唨捨入誤差嘅後果。喺歐拉法嘅𨂾數n入便,捨入誤差大致係εyn幅值,其中ε係嘅機器ε。假設捨入誤差大致都係相同大細嘅,喺N𨂾之後佮埋嘅捨入誤差係大致Nεy0若果所有錯誤指嗮相同方向。由於𨂾數戥𨂾長h成反比,所以總捨入誤差戥ε/ h成正比。但係,實質上,所有捨入誤差唔可能都指嗮同一方向。相反,若果假設捨入誤差係獨立嘅隨機變數,噉預期嘅總捨入誤差正比於[19]

所以,對於極細嘅𨂾長值,截尾誤差會好細,但係捨入誤差嘅影響可能會好大。若果喺Euler法嘅公式裏便使埋補償求和,可以好容易噉避免捨入誤差嘅大部分影響。 [20]

修改同擴展

[編輯]

後褪Euler法係對Euler法嘅簡單修改,消除唨上一部分入便提到嘅穩定性問題:

噉樣戥(標准或順向嘅)歐拉法唔同之處係個函數喺一𨂾嘅終點得到評估,而唔係起點。後褪Euler法係一個隱式方法(所以又喊做隱式Euler法),代表後褪Euler法嘅公式喺兩便都有,所以喺應用後褪歐拉法嗰陣,必須求解一個方程。噉樣令到實施更加昂貴。

其他啲對Euler法有助於穩定性嘅修改,產生唨指數Euler法半隱式Euler法

更複雜嘅方法可以實現更高嘅階數(以及更高嘅準確性)。一種可能性係使更多函數評估值(function evaluations)。對此,用返本文提過嘅中點法進行說明:

噉樣導向一系列Runge-Kutta法

另一種可能性係使多尐過去值,似兩步Adams–Bashforth方法所示:

噉樣導向一系列線性多步法。仲有其他一啲修改,啲修改使用壓縮感測(compressive sensing)嘅技巧來最大程度噉減少內存使用。[21]

流行文化

[編輯]

喺電影《NASA無名英雌》入便,Katherine Goble借助歐拉法來計算宇航員John Glenn由地球軌道嘅再入[22]

睇埋

[編輯]

[編輯]
  1. Butcher 2003, p. 45; Hairer, Nørsett & Wanner 1993, p. 35
  2. Atkinson 1989, p. 342; Butcher 2003, p. 60
  3. Butcher 2003, p. 45; Hairer, Nørsett & Wanner 1993, p. 36
  4. Butcher 2003, p. 3; Hairer, Nørsett & Wanner 1993, p. 2
  5. See also Atkinson 1989, p. 344
  6. Hairer, Nørsett & Wanner 1993, p. 40
  7. Atkinson 1989, p. 342; Hairer, Nørsett & Wanner 1993, p. 36
  8. Atkinson 1989, p. 342
  9. Atkinson 1989, p. 343
  10. Butcher 2003, p. 60
  11. Atkinson 1989, p. 342
  12. Stoer & Bulirsch 2002, p. 474
  13. Atkinson 1989, p. 344
  14. Butcher 2003, p. 49
  15. Atkinson 1989, p. 346; Lakoba 2012
  16. Iserles 1996, p. 7
  17. Butcher 2003, p. 63
  18. Butcher 2003, p. 70; Iserles 1996, p. 57
  19. Butcher 2003
  20. Butcher 2003
  21. Unni, M. P.; Chandra, M. G.; Kumar, A. A. (March 2017). "Memory reduction for numerical solution of differential equations using compressive sensing". 2017 IEEE 13th International Colloquium on Signal Processing Its Applications (CSPA): 79–84. doi:10.1109/CSPA.2017.8064928. ISBN 978-1-5090-1184-1.
  22. Khan, Amina. "Meet the 'Hidden Figures' mathematician who helped send Americans into space". Los Angeles Times. 喺12 February 2017搵到.

[編輯]

連出去

[編輯]