歐拉法
喺數學同計算科學入便,歐拉法(亦都叫做正向歐拉法)係一種一階數值過程,攞嚟求解啲具有畀定初始值嘅常微分方程(ODE)嘅。歐拉法係對常微分方程進行數值積分嘅最基本嘅顯式方法,亦都係最簡單嘅Runge-Kutta法。歐拉法以歐拉個名命名,佢喺佢本《計算積分學》(1768-1870年出版)一書入便有處理到。 [1]
歐拉法係一階方法,代表佢個局部誤差(每𨂾誤差)戥𨂾長嘅平方成正比、全局誤差(畀定時間嘅誤差)戥𨂾長成正比。歐拉法通常作為構造啲更複雜方法嘅基礎,譬如預測器-校正器法。
非正式嘅幾何描述
編輯考慮計算某條未知曲線個形狀;佢從畀定點開始、滿足一條畀定微分方程。基於呢層,微分方程即係一種公式,一旦計出任意一粒點個位置,透過呢條公式就計得出條曲線上嘅切線斜率嘅。
個惗法係,雖然條曲線初初都係未知嘅,但佢個起點 係已知嘅(請參見右上方幅圖)。然之後,從微分方程就計得出條曲線喺 度嘅斜率同埋切線。
沿住切線向上延伸多一點 。沿住呢一細𨂾,個斜率唔會變化太大,所以 好接近條曲線。若果假設個 仲喺曲線上,戥上述喺點 相同嘅推理過程都繼續使得。幾𨂾之後,就計出多邊形曲線 。好多時,呢條曲線戥原始未知曲線嘅走差唔會太大;若果𨂾長夠細而且計算間隔有限,噉可以令到兩條曲線之間嘅誤差細啲: [2]
揀一個值 對於每一𨂾同埋設置 。噉時,從 到 係: [3]
值 係ODE解嘅近似值,喺 : 嗰陣。 Euler法係顯式嘅,即對於 ,解 係 嘅顯式函數。
就算Euler法集成唨一階ODE,但係任何N階嘅ODE都可以表示為啲一階ODE嘅系統。
爲唨處理以下方程:
- ,
可以引入輔助變量 並得到等式:
呢個係變量 入便嘅一階系統,可以用Euler法處理,或者實質上可以用任一其他啲都用於一階系統嘅方案處理。 [4]
例子
編輯畀定初始值問題
諗住用歐拉法來逼近 。 [5]
令𨂾長等於1(h = 1)
編輯歐拉法係
所以首先必須計算 。喺條簡單嘅微分方程入便,函數 由定義 。有初始值:
透過執行以上步驟,就搵到條切線嘅斜率、喺點 戥方程解條曲線相切嘅 。考慮返斜率定義為 變化量除以 變化量 ,或者 ,
下一步係將上述值乘上𨂾長 ,呢度等於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
使埋其他𨂾長
編輯好似引言入便所述,若果𨂾長 細啲,噉歐拉法亦都準確啲。下表顯示唨唔同𨂾長嘅結果。第一行對應於上一節入便個示例,第二行如圖所示。
𨂾長 | 歐拉法結果 | 誤差 |
---|---|---|
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]相似嘅計算推得出中點法同埋後褪歐拉法。
最後,可以積分條微分方程式由 到 ,並應用微積分基本定理得到:
再透過左手矩形法(只用到一個矩形)來近似積分:
局部截尾誤差
編輯歐拉法嘅局部截尾誤差(LTE)係行親一𨂾造成嘅誤差,即係一𨂾之後數值解 戥 嗰陣嘅確切解之間嘅差。數值解由下式畀出
對於確切嘅解,使返上便嘅#推導部分入便提到嘅泰勒展開式:
歐拉法引入嘅局部截尾誤差由以下等式之間嘅差異畀出:
呢個結果有效,若果 有有界三階導數。[10]
噉樣表明對於細 ,局部截尾誤差正比於 ,即會令到Euler法對於細啲嘅 準確性降低,而唔似Runge-Kutta法同埋線性多步法之類其他啲高階竅,因為後者啲局部截尾誤差戥𨂾長嘅高啲嘅冪成正比。
一條略有唔同嘅局部截尾誤差公式,可以透過對於泰勒定理入便啲其餘項、用返拉格朗日形式獲得。若果 有一個連續嘅二階導數,噉就存在一個 令到:
喺上便嘅誤差表達式裏便,未知精確解 嘅二階導數可以用一條表達式、涉及微分方程右側嘅嚟代替。根據 方程有:
全局截尾誤差
編輯全局截尾誤差(GTE)係固定時間嘅誤差 ,但係係經過好多𨂾步後,先從初始時間嗰度開始計。全局截尾誤差係每一𨂾入便整出嘅嗰啲局部截尾誤差嘅累積影響。 [13]𨂾數好容易確定為 ,之正比於 ,而且每一𨂾裏便整出嘅誤差正比於 (請參閱上一節)。所以,預期當中,全局截尾誤差會正比於 。 [14]
呢種直觀嘅推理可以精確化。若果個解 有有界二階導數、若果 喺佢入便個第二個參數係Lipschitz連續嘅,噉個全局截尾誤差由條式得出:
其中 係喺畀定間隔下、 個二階導數嘅上限,而 係 嘅Lipschitz常數 。 [15]
呢個界限嘅精確形式幾乎冇乜實際意義,因為喺大多數情況之下,界限遠遠高估唨Euler法造成嘅實際誤差。 [16]重要嘅衹係佢表明到個全局截尾誤差(大約)正比於 。所以,歐拉法着認為係一階嘅。 [17]
數值穩定性
編輯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]
睇埋
編輯- Crank–Nicolson法(Crank–Nicolson method)
- 梯度下降同樣使用有限𨂾長,搵到函數嘅最細值
- Runge-Kutta法同埋個列表
- 線性多步法(Linear multistep method)
- 數值積分(攞嚟計算定積分)
- 常微分方程嘅數值方法(Numerical methods for ordinary differential equations)
註
編輯- ↑ Butcher 2003, p. 45; Hairer, Nørsett & Wanner 1993, p. 35
- ↑ Atkinson 1989, p. 342; Butcher 2003, p. 60
- ↑ Butcher 2003, p. 45; Hairer, Nørsett & Wanner 1993, p. 36
- ↑ Butcher 2003, p. 3; Hairer, Nørsett & Wanner 1993, p. 2
- ↑ See also Atkinson 1989, p. 344
- ↑ Hairer, Nørsett & Wanner 1993, p. 40
- ↑ Atkinson 1989, p. 342; Hairer, Nørsett & Wanner 1993, p. 36
- ↑ Atkinson 1989, p. 342
- ↑ Atkinson 1989, p. 343
- ↑ Butcher 2003, p. 60
- ↑ Atkinson 1989, p. 342
- ↑ Stoer & Bulirsch 2002, p. 474
- ↑ Atkinson 1989, p. 344
- ↑ Butcher 2003, p. 49
- ↑ Atkinson 1989, p. 346; Lakoba 2012
- ↑ Iserles 1996, p. 7
- ↑ Butcher 2003, p. 63
- ↑ Butcher 2003, p. 70; Iserles 1996, p. 57
- ↑ Butcher 2003
- ↑ Butcher 2003
- ↑ 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.
- ↑ Khan, Amina. "Meet the 'Hidden Figures' mathematician who helped send Americans into space". Los Angeles Times. 喺12 February 2017搵到.
攷
編輯- Atkinson, Kendall A. (1989). An Introduction to Numerical Analysis (第2版). New York: John Wiley & Sons. ISBN 978-0-471-50023-0.
- Ascher, Uri M.; Petzold, Linda R. (1998). Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations. Philadelphia: Society for Industrial and Applied Mathematics. ISBN 978-0-89871-412-8.
- Butcher, John C. (2003). Numerical Methods for Ordinary Differential Equations. New York: John Wiley & Sons. ISBN 978-0-471-96758-3.
- Hairer, Ernst; Nørsett, Syvert Paul; Wanner, Gerhard (1993). Solving ordinary differential equations I: Nonstiff problems. Berlin, New York: Springer-Verlag. ISBN 978-3-540-56670-0.
- Iserles, Arieh (1996). A First Course in the Numerical Analysis of Differential Equations. Cambridge University Press. ISBN 978-0-521-55655-2.
- Stoer, Josef; Bulirsch, Roland (2002). Introduction to Numerical Analysis (第3版). Berlin, New York: Springer-Verlag. ISBN 978-0-387-95452-3.
- Lakoba, Taras I. (2012), Simple Euler method and its modifications (PDF) (Lecture notes for MATH334), University of Vermont, 喺29 February 2012搵到
- Unni, M P. (2017). "Memory reduction for numerical solution of differential equations using compressive sensing". 2017 IEEE 13th International Colloquium on Signal Processing & its Applications (CSPA). IEEE CSPA. pp. 79–84. doi:10.1109/CSPA.2017.8064928. ISBN 978-1-5090-1184-1.
連出去
編輯- 維基同享中有關歐拉法嘅資源類</img>
- Rosetta Code用唔同語言實現嘅Euler法實現
- Hazewinkel, Michiel, 編 (2001), "Euler method", Encyclopaedia of Mathematics, Kluwer Academic Publishers, ISBN 978-1556080104