歐拉法

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

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

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

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

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

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

揀一個值h對於每一𨂾同埋設置t_n = t_0 + nh 。噉時，從t_n到t_{n+1} = t_n + h係： 

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

值y_n係ODE解嘅近似值，喺t_n ： y_n \approx y(t_n) 嗰陣。 Euler法係顯式嘅，即對於i \leq n，解y_{n+1}係y_i嘅顯式函數。

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

爲唨處理以下方程：

:  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)入便嘅一階系統，可以用Euler法處理，或者實質上可以用任一其他啲都用於一階系統嘅方案處理。 

## 例子
畀定初始值問題

: y'=y, \quad y(0)=1, 

諗住用歐拉法來逼近y(4) 。 

### 令𨂾長等於1（''h'' = 1）

歐拉法係

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

所以首先必須計算f(t_0, y_0) 。喺條簡單嘅微分方程入便，函數f由定義f(t,y) = y 。有初始值：

:  f(t_0,y_0) = f(0,1) = 1. \qquad \qquad

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

下一步係將上述值乘上𨂾長h ，呢度等於1：

:  h \cdot f(y_0) = 1 \cdot 1 = 1. \qquad \qquad

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

:  y_0 + hf(y_0) = y_1 = 1 + 1 \cdot 1 = 2. \qquad \qquad

重複上述步驟即可搵到 y_2、y_3 、y_4 。

:  \begin{align}
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_4 &= y_3 + hf(y_3) = 8 + 1 \cdot 8 = 16.
\end{align} 

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

個計算嘅結論係 y_4 = 16  。微分方程嘅精確解係 y(t) = e^t  ， 所以 y(4) = e^4 \approx 54.598  。就算喺呢種𨂾長值h大嘅特定情況之下，歐拉法嘅逼近唔點精確，之不過如圖所示，演算法個表現喺定性上係啱嘅。

### MATLAB代碼示例

```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編程語言編寫嘅示例代碼。
```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細啲，噉歐拉法亦都準確啲。下表顯示唨唔同𨂾長嘅結果。第一行對應於上一節入便個示例，第二行如圖所示。

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

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

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

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

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

:  y(t_0 + h) = y(t_0) + h y'(t_0) + \frac{1}{2}h^2 y*(t_0) + O(h^3). </math>

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

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

: <math> y'(t_0) \approx \frac{y(t_0+h) - y(t_0)}{h} </math>

喺微分方程入便<math>y' = f(t,y)</math> 。同樣，噉樣都推得出歐拉法。 相似嘅計算推得出[[中點法]]同埋[[後褪歐拉法]]。

最後，可以積分條微分方程式由<math>t_0</math>到<math>t_0 + h</math>，並應用[[微積分基本定理|微積分]]基本定理得到：

: <math> y(t_0+h) - y(t_0) = \int_{t_0}^{t_0+h} f(t,y(t)) \,\mathrm{d}t. </math>

再透過左手[[矩形法]]（只用到一個矩形）來近似積分：

: <math> \int_{t_0}^{t_0+h} f(t,y(t)) \,\mathrm{d}t \approx h f(t_0, y(t_0)). </math>

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

== 局部截尾誤差 ==
歐拉法嘅[[局部截尾誤差]]（LTE）係行親一𨂾造成嘅誤差，即係一𨂾之後數值解<math>y_1</math>戥<math>t_1 = t_0+h</math>嗰陣嘅確切解之間嘅差。數值解由下式畀出

: <math> y_1 = y_0 + h f(t_0, y_0). \quad</math>

對於確切嘅解，使返上便嘅[[#推導]]部分入便提到嘅泰勒展開式：

: <math> 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). </math>

呢個結果有效，若果<math>y</math>有有界三階導數。

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

一條略有唔同嘅局部截尾誤差公式，可以透過對於泰勒定理入便啲其餘項、用返拉格朗日形式獲得。若果<math>y</math>有一個連續嘅二階導數，噉就存在一個<math>\xi \in [t_0,t_0+h]</math>令到：

: <math> \mathrm{LTE} = y(t_0 + h) - y_1 = \frac{1}{2} h^2 y*(\xi).  

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

: y*(t_0) = {\partial f\over\partial t}(t_0, y(t_0)) + {\partial f\over\partial y}(t_0, y(t_0)) \, f(t_0, y(t_0)).</math> 

== 全局截尾誤差 ==
[[全局截尾誤差]]（GTE）係固定時間嘅誤差<math>t</math> ，但係係經過好多𨂾步後，先從初始時間嗰度開始計。全局截尾誤差係每一𨂾入便整出嘅嗰啲局部截尾誤差嘅累積影響。 𨂾數好容易確定為<math>(t-t_0)/h</math> ，之正比於<math>1/h</math> ，而且每一𨂾裏便整出嘅誤差正比於<math>h^2</math> （請參閱上一節）。所以，預期當中，全局截尾誤差會正比於<math>h</math> 。 

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

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

其中<math>M</math>係喺畀定間隔下、<math>y</math>個二階導數嘅上限，而<math>L</math>係<math>f</math>嘅Lipschitz常數 。 

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

== 數值穩定性 ==

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

: <math> y' = -2.3y, \qquad y(0) = 1. </math>

確切解係<math>y(t) = e^{-2.3t}</math> ，隨住<math>t \to \infty</math> 而到零。但係，若果歐拉法以𨂾長<math>h=1</math>應用於呢條式，噉個數值解喺定性上都係錯嘅：個數值解會振盪並增長（見圖），即唔穩定噉解。又如，若果𨂾長細啲、<math>h = 0.7</math> ，噉個數值解又確實會衰減為零。

若果將歐拉法應用於線性方程<math>y' = k y</math> ，噉呢個乘積嘅數值解會唔穩定、若果<math>hk</math>乘積喺以下區域之外：

: <math> \{ z \in \mathbf{C} \mid |z+1| \le 1 \}, </math>

似右圖所示，嗰柵區域叫做（線性）*'穩定性區域*'。 喺呢個例子入便<math>k</math>係-2.3，所以若果<math>h = 1</math>然之後<math>hk = -2.3</math>，喺穩定區域之外，所以數值解唔穩定。

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

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

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

## 修改同擴展
後褪Euler法係對Euler法嘅簡單修改，消除唨上一部分入便提到嘅穩定性問題：

:  y_{n+1} = y_n + h f(t_{n+1}, y_{n+1}). 

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

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

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

:  y_{n+1} = y_n + h f \Big( t_n + \tfrac12 h, y_n + \tfrac12 h f(t_n, y_n) \Big) 。

噉樣導向一系列Runge-Kutta法。

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

:  y_{n+1} = y_n + \tfrac32 h f(t_{n}, y_{n}) - \tfrac12 h f(t_{n-1}, y_{n-1}). 

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

## 流行文化
喺電影《NASA無名英雌》入便，Katherine Goble借助歐拉法來計算宇航員John Glenn由地球軌道嘅再入。 

## 睇埋

* Crank–Nicolson法（Crank–Nicolson method）
* 梯度下降同樣使用有限𨂾長，搵到函數嘅最細值
* Runge-Kutta法同埋
* 線性多步法（Linear multistep method）
* 數值積分（攞嚟計算定積分）
* 常微分方程嘅數值方法（Numerical methods for ordinary differential equations）

## 註

## 攷

* 
* 
* 
* 
* 
* 
* 
* 

## 連出去

* 
* Rosetta Code用唔同語言實現嘅Euler法實現
*