✨Phương pháp Euler

thumb|Minh họa phương pháp Euler. Đường cong chưa biết có màu xanh da trời và lời giải gần đúng của nó là đường nhiều cạnh màu đỏ. Trong toán học và khoa học máy tính, phương pháp Euler là một phương pháp số bậc một để giải các phương trình vi phân thường (ODEs) với giá trị ban đầu cho trước. Nó là phương pháp hiện (explicit) cơ bản nhất cho việc tính tích phân số của các phương trình vi phân thường và là phương pháp Runge-Kutta đơn giản nhất. Phương pháp Euler được đặt theo tên của Leonhard Euler, người đã đề cập đến phương pháp này trong cuốn sách Institutionum calculi integralis của ông (xuất bản 1768-1770).

Phương pháp Euler là một phương pháp bậc một, có nghĩa là sai số cục bộ (sai số mỗi bước) tỷ lệ thuận với bình phương của kích thước bước, và sai số tổng thể (sai số tại một thời điểm nào đó) tỷ lệ thuận với kích thước bước. Phương pháp Euler thường phục vụ như là cơ sở để xây dựng các phương pháp phức tạp hơn, ví dụ như, phương pháp Dự đoán- Hiệu chỉnh.

Mô tả hình học phi chính thức

Xem xét bài toán tính toán hình dạng của một đường cong chưa biết bắt đầu tại một điểm cho trước và thỏa mãn một phương trình vi phân nào đó. Ở đây, phương trình vi phân có thể được coi như là một công thức mà nhờ nó độ dốc của đường tiếp tuyến với đường cong có thể được tính toán tại điểm bất kỳ nào trên đường cong, một khi vị trí của điểm đó đã được tính toán.

Ý tưởng là trong khi đường cong ban đầu chưa biết, điểm khởi đầu của nó được biểu thị bởi $A\_0,$ là đã biết (xem hình phía trên bên phải). Thì, từ phương trình vi phân, độ dốc của đường cong tại $A\_0$ có thể được tính toán, và vì vậy, có thể tìm được đường tiếp tuyến.

Đi một bước nhỏ dọc theo đường tiếp tuyến đến một điểm $A\_1.$ Dọc bước nhỏ này, độ dốc không thay đổi quá nhiều, vì vậy $A\_1$ sẽ gần với đường cong. Nếu chúng ta coi $A\_1$ vẫn còn nằm trên đường cong, cùng một lý luận như đối với điểm $A\_0$ có thể tính được độ dốc và tiếp tuyến của đường cong tại $A\_1$. Sau vài bước, một đường cong đa giác $A\_0A\_1A\_2A\_3\backslash dots$ được tìm ra. Nói chung, đường cong này không phân kỳ quá xa khỏi đường cong chưa biết ban đầu, và sai số giữa hai đường cong có thể được làm nhỏ nếu kích thước bước là đủ nhỏ và khoảng thời gian tính toán là hữu hạn. :$y\text{'}(t)\; =\; f(t,y(t)),\; \backslash qquad\; \backslash qquad\; y(t\_0)=y\_0.$

Chọn một giá trị $h$ cho kích thước của mỗi bước và đặt $t\_n\; =\; t\_0\; +\; nh$. Bây giờ, một bước của phương pháp Euler từ $t$n tới $t${n+1} = t_n + h là

:$y\_\{n+1\}\; =\; y\_n\; +\; hf(t\_n,y\_n).$

Giá trị của $y\_n$ là một lời giải gần đúng của phương trình ODE tại thời điểm $t\_n$: $y\_n\; \backslash approx\; y(t$n). Phương pháp Euler là phương pháp hiện, nghĩa là lời giải $y${n+1} là một hàm hiện của $y\_i$ với $i\; \backslash leq\; n$.

Trong khi phương pháp Euler tích phân một ODE bậc nhất, ODE bất kỳ bậc N có thể được biểu diễn như là một ODE bậc nhất: để xử lý phương trình

:$y^\{(N)\}(t)\; =\; f(t,\; y(t),\; y\text{'}(t),\; \backslash ldots,\; y^\{(N-1)\}(t))$,

chúng ta giới thiệu các biến phụ $z\_1(t)=y(t),\; z\_2(t)=y\text{'}(t),\backslash ldots,\; z\_N(t)=y^\{(N-1)\}(t)$ và có được phương trình tương đương

:$\backslash mathbf\{z\}\text{'}(t)\; =\; \backslash 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}

Đây là một hệ bậc nhất với biến $\backslash mathbf\{z\}(t)$ và có thể được giải bằng phương pháp Euler hoặc bằng bất kỳ lược đồ nào được sử dụng để giải các hệ bậc nhất.

Ví dụ

Cho bài toán giá trị ban đầu

:$y\text{'}=y,\; \backslash quad\; y(0)=1,$

chúng ta muốn sử dụng phương pháp Euler để xấp xỉ $y(4)$.

Sử dụng kích thước bước bằng 1 (h = 1)

thumb|Minh họa tích phân số cho phương trình $y\text{'}=y,\; y(0)=1.$ Màu xanh da trời là phương pháp Euler; màu xanh lá cây là[[midpoint method|phương pháp điểm giữa; màu đỏ là lời giải chính xác, $y=e^t.$ Kích thước bước là h = 1.0.]] Phương pháp Euler là:

:$y\_\{n+1\}\; =\; y\_n\; +\; hf(t\_n,y\_n).\; \backslash qquad\; \backslash qquad$

vì vậy trước tiên chúng ta phải tính toán $f(t\_0,\; y\_0)$. Trong phương trình vi phân đơn giản này, hàm $f$ được định nghĩa bởi $f(t,y)\; =\; y$. Chúng ta có

:$f(t\_0,y\_0)\; =\; f(0,1)\; =\; 1.\; \backslash qquad\; \backslash qquad$

Bằng cách thực hiện bước trên, chúng ta đã tìm được độ dốc của đường thẳng tiếp tuyến với đường cong lời giải tại điểm $(0,1)$.

Nhớ lại rằng độ dốc được định nghĩa là sự thay đổi trong $y$ chia cho sự thay đổi trong $t$, tức là $\backslash Delta\; y/\; \backslash Delta\; t$.

Bước tiếp theo là nhân giá trị trên với kích thước bước $h$, cái chúng ta lấy bằng một ở đây:  :$h\; \backslash cdot\; f(y\_0)\; =\; 1\; \backslash cdot\; 1\; =\; 1.\; \backslash qquad\; \backslash qquad$

Bởi vì kích thước bước là sự thay đổi trong $t$, khi chúng ta nhân kích thước bước và độ dốc của tiếp tuyến, chúng ta có được sự thay đổi trong giá trị $y$. Giá trị này sau đó sẽ được thêm vào giá trị $y$ ban đầu để có được giá trị tiếp theo được sử dụng để tính toán.

:$y\_0\; +\; hf(y\_0)\; =\; y\_1\; =\; 1\; +\; 1\; \backslash cdot\; 1\; =\; 2.\; \backslash qquad\; \backslash qquad$

Các bước trên sẽ được lặp đi lặp lại để tìm ra $y\_2$, $y\_3$ và $y\_4$.

:

Do tính chất lặp đi lặp lại của thuật toán này, nên tổ chức các bước tính toán dưới dạng biểu đồ, như thấy bên dưới, để tránh tạo ra các sai sót.

:

Đáp án cuối cùng là $y\_4\; =\; 16$. Lời giải chính xác của phương trình vi phân là $y(t)\; =\; e^t$, vì vậy $y(4)\; =\; e^4\; \backslash approx\; 54.598$. Như vậy, lời giải gần đúng của phương pháp Euler không thật sự tốt trong trường hợp này. Tuy nhiên, như hình vẽ cho thấy, ứng xử của nó là đúng về mặt định tính.

Sử dụng các kích thước bước khác

thumb|Với cùng minh họa cho trường hợp h = 0.25. Như đã nói trong phần giới thiệu, phương pháp Euler sẽ chính xác hơn nếu kích thước bước $h$ nhỏ hơn. Bảng dưới đây cho thấy kết quả với các kích thước bước khác nhau. Hàng trên cùng tương ứng với ví dụ trong phần trước, và hàng thứ hai được minh họa trong hình.

:

Sai số được ghi ở cột cuối cùng của bảng là sự khác biệt giữa lời giải chính xác tại $t\; =\; 4$ và lời giải gần đúng Euler. Ở dưới cùng của bảng, kích thước bước bằng một nửa kích thước bước ở hàng trên nó, và sai số cũng bằng khoảng một nửa sai số trong dòng trên nó. Điều này cho thấy rằng sai số gần như là tỷ lệ với kích thước bước, ít nhất là đối với các giá trị kích thước bước tương đối nhỏ. Điều này nói chung, cũng đúng cho các phương trình khác; xem phần Sai số cắt cụt tổng thể để biết thêm chi tiết.

Các phương pháp khác, chẳng hạn như phương pháp điểm giữa cũng được minh họa trên các hình, và có vẻ chính xác hơn: sai số của phương pháp điểm giữa gần như tỷ lệ thuận với bình phương của kích thước bước. Vì lý do này, phương pháp Euler được gọi là một phương pháp bậc nhất, trong khi phương pháp điểm giữa là phương pháp bậc hai.

Chúng ta có thể ngoại suy từ bảng trên rằng kích thước bước cần thiết để có được một đáp án chính xác đến ba chữ số thập phân là khoảng 0,00001, có nghĩa là chúng ta cần 400.000 bước. Con số bước rất lớn này đòi hỏi một chi phí tính toán cao. Vì lý do này, người ta thường sử dụng các phương pháp bậc cao, thay thế cho phương pháp Euler, như là các phương pháp Runge-Kutta hoặc các phương pháp đa bước tuyến tính, đặc biệt là nếu muốn đạt được độ chính xác cao.

Nguồn gốc

Phương pháp Euler có thể được rút ra theo một số cách. Cách thứ nhất, đó là mô tả hình học đã được đề cập ở trên.

Một khả năng khác là xem xét mở rộng Taylor của hàm $y$ around $t\_0$:  :$y(t\_0\; +\; h)\; =\; y(t\_0)\; +\; h\; y\text{'}(t\_0)\; +\; \backslash frac\{1\}\{2\}h^2\; y\text{'}\text{'}(t\_0)\; +\; O(h^3).$

Phương trình vi phân $y\text{'}=f(t,y)$. Nếu phương trình vi phân này thay vào trong mở rộng Taylor và các số hạng bậc hai và bậc cao hơn bị bỏ qua, thì sẽ có được phương pháp Euler. Mở rộng Taylor được sử dụng dưới đây để phân tích sai số của phương pháp Euler, và nó có thể được mở rộng để đạt được các phương pháp Runge-Kutta.

Một cách khác là để thay thế công thức sai phân hữu hạn tiếp tới cho đạo hàm, :$y\text{'}(t\_0)\; \backslash approx\; \backslash frac\{y(t\_0+h)\; -\; y(t\_0)\}\{h\}$

trong phương trình vi phân $y\text{'}\; =\; f(t,y)$. Một lần nữa, điều này đem lại phương pháp Euler.^[8] Theo cách tương tự sẽ dẫn đến quy tắc điểm giữa và phương pháp Euler lùi lại.

Cuối cùng, ta có thể tích phân phương trình vi phân từ $t\_0$ tới $t\_0\; +\; h$ và áp dụng các định lý cơ bản của tích phân và vi phân để có được:

:$y(t\_0+h)\; -\; y(t$0) = \int{t_0}^{t_0+h} f(t,y(t)) \,\mathrm{d}t.

Bây giờ xấp xỉ tích phân bằng phương pháp hình chữ nhật bên trái:

:$\backslash int\_\{t\_0\}^\{t\_0+h\}\; f(t,y(t))\; \backslash ,\backslash mathrm\{d\}t\; \backslash approx\; h\; f(t\_0,\; y(t\_0)).$

Kết hợp cả hai phương trình, ta lại tìm thấy phương pháp Euler một lần nữa. Cách tiếp cận này có thể được tiếp tục để đi đến nhiều phương pháp đa bước tuyến tính khác.

Sai số cắt cụt cục bộ

Sai số cắt cụt cục bộ của phương pháp Euler là sai số trong một bước duy nhất. Đó là sự khác biệt giữa lời giải số sau một bước, $y\_1$, và lời giải chính xác tại thời điểm $t\_1\; =\; t\_0+h$. Lời giải số được cho bởi

:$y\_1\; =\; y\_0\; +\; h\; f(t\_0,\; y\_0).\; \backslash quad$

Đối với lời giải chính xác, chúng ta sử dụng mở rộng Taylor được đề cập trong phần Nguồn gốc phía trên:

:$y(t\_0\; +\; h)\; =\; y(t\_0)\; +\; h\; y\text{'}(t\_0)\; +\; \backslash frac\{1\}\{2\}h^2\; y\text{'}\text{'}(t\_0)\; +\; O(h^3).$

Sai số cắt cụt cục bộ (LTE) của phương pháp Euler được cho bởi sự khác biệt giữa các phương trình này:

:$\backslash mathrm\{LTE\}\; =\; y(t\_0\; +\; h)\; -\; y\_1\; =\; \backslash frac\{1\}\{2\}\; h^2\; y\text{'}\text{'}(t\_0)\; +\; O(h^3).$

Kết quả này là hợp lý nếu $y$ có một đạo hàm bậc ba bị chặn (bounded).

Điều này cho thấy rằng đối với $h$ nhỏ, các sai số cắt cụt cục bộ xấp xỉ tỷ lệ thuận với $h^2$ do đó làm cho phương pháp Euler kém chính xác (đối với $h$ nhỏ) hơn so với các phương pháp bậc cao khác như các phương pháp Runge-Kutta và các phương pháp đa bước tuyến tính, mà sai số cắt cụt cục bộ tỷ lệ thuận với một số mũ cao hơn của kích thước bước.

Một cách xây dựng công thức hơi khác cho sai số cắt cụt cục bộ là sử dụng dạng thức Lagrange cho số hạng còn lại trong định lý Taylor. Nếu $y$ có đạo hàm bậc hai liên tục, thì tồn tại một $\backslash xi\; \backslash in\; [t\_0,t\_0+h]$ mà :$\backslash mathrm\{LTE\}\; =\; y(t\_0\; +\; h)\; -\; y\_1\; =\; \backslash frac\{1\}\{2\}\; h^2\; y\text{'}\text{'}(\backslash xi).$

Trong các biểu thức sai số trên, đạo hàm bậc hai của lời giải chính xác chưa biết $y$ có thể được thay thế bằng một biểu thức ở phía bên phải của phương trình vi phân. Thật vậy, từ phương trình $y\text{'}=f(t,y)$ ta có

:$y\text{'}\text{'}(t\_0)\; =\; \{\backslash partial\; f\backslash over\backslash partial\; t\}(t\_0,\; y(t\_0))\; +\; \{\backslash partial\; f\backslash over\backslash partial\; y\}(t\_0,\; y(t\_0))\; \backslash ,\; f(t\_0,\; y(t\_0)).$

Sai số cắt cụt tổng thể

Sai số cắt cụt tổng thể là sai số tại một thời điểm cố định $t$, sau nhiều bước nhiều phương pháp cần phải thực hiện để đạt được thời điểm đó từ thời điểm ban đầu. Sai số cắt cụt tổng thể là tích lũy của các sai số cắt cụt cục bộ đã phạm phải trong mỗi bước trước đó. Số lượng các bước được dễ dàng xác định là $(t-t\_0)/h$, tỷ lệ thuận với $1/h$, và sai số đã phạm phải trong mỗi bước tỷ lệ thuận với $h^2$(xem phần trước). Vì vậy, mong đợi rằng sai số cắt cụt tổng thể sẽ tỷ lệ thuận với $h$.

Lý luận trực quan này có thể được chứng minh là chính xác. Nếu lời giải $y$ có đạo hàm bậc hai bị chặn (bounded) và $f$ là Lipschitz liên tục trong đối số thứ hai của nó, thì sai số cắt cụt tổng thể (GTE) được bao (bounded) bởi: :$|\backslash text\{GTE\}|\; \backslash le\; \backslash frac\{hM\}\{2L\}(e^\{L(t-t\_0)\}-1)\; \backslash qquad\; \backslash qquad$ trong đó $M$ là một giới hạn trên cho đạo hàm bậc hai của $y$ trên khoảng thời gian nhất định nào đó và $L$ là hằng số Lipschitz của $f$.

Dạng thức chính xác của giới hạn này ít quan trọng trong thực tế, trong hầu hết các trường hợp, giới hạn này quá lớn so với sai số thực sự phạm phải bởi phương pháp Euler.^[16] Điều quan trọng là nó cho thấy rằng sai số cắt cụt tổng thể (một cách gần đúng) tỷ lệ thuận với. Vì lý do này, phương pháp Euler được cho/ gọi là bậc nhất.

Sự ổn định số của phương pháp Euler

thumb|Lời giải của phương trình $y\text{'}\; =\; -2.3y$ được tính toán bằng phương pháp Euler với kích thước bước $h=1$ (các hình vuông màu xanh da trời) và $h=0.7$ (các hình tròn màu đỏ). Đường cong màu đen biểu diễn lời giải chính xác. Phương pháp Euler cũng có thể không ổn định về mặt phương pháp số, đặc biệt là đối với các phương trình cứng, có nghĩa là lời giải số tăng rất nhanh trong khi lời giải chính xác không (tăng). Điều này có thể được minh họa bằng cách sử dụng phương trình tuyến tính :$y\text{'}\; =\; -2.3y,\; \backslash qquad\; y(0)\; =\; 1.$ thumb|Hình tròn lớn màu hồng biểu diễn vùng ổn định đối với phương pháp Euler Lời giải chính xác là $y(t)\; =\; e^\{-2.3t\}$, phân rã về không khi $t\; \backslash to\; \backslash infty$. Tuy nhiên, nếu phương pháp Euler được áp dụng cho phương trình này với kích thước bước $h=1$, thì lời giải số là sai về mặt định tính: nó dao động và tăng (xem hình). Đây là những gì có nghĩa là không ổn định. Nếu kích thước bước nhỏ hơn được sử dụng, ví dụ $h=0.7$, thì lời giải số không phân rã về không.

Nếu phương pháp Euler được áp dụng cho phương trình tuyến tính $y\text{'}\; =\; k\; y$, thì lời giải số không ổn định nếu tích số $hk$ nằm bên ngoài vùng :$\{\; z\; \backslash in\; \backslash mathbf\{C\}\; \backslash mid\; |z+1|\; \backslash le\; 1\; \},$ được minh họa ở phía bên phải. Vùng này được gọi là vùng không ổn định (tuyến tính). Trong ví dụ, $k$ bằng -2,3, vì vậy nếu $h=1$ thì $hk\; =\; -2,3$ tức là nằm bên ngoài vùng ổn định, và do đó lời giải số là không ổn định.

Sự hạn chế này —cùng với việc hội tụ sai số chậm với h— làm cho phương pháp Euler không được sử dụng thường xuyên, ngoại trừ như một ví dụ đơn giản của tích phân số.

Các sai số làm tròn

Thảo luận từ trên đến giờ đã bỏ qua những hậu quả của sai số làm tròn. Trong bước n của phương pháp Euler, sai số làm tròn là xấp xỉ độ lớn εy_n trong đó ε là Machine epsilon (giới hạn trên của sai số tương đối do làm tròn trong số học điểm nổi). Giả sử rằng các sai số làm tròn tất cả có kích thước xấp xỉ như nhau, sai số làm tròn tổng hợp trong N bước là xấp xỉ _N_εy₀ nếu tất cả các sai số chỉ về cùng hướng. Bởi vì số lượng bước tỉ lệ nghịch với kích thước bước h, tổng sai số làm tròn tỷ lệ thuận với ε / h. Trong thực tế, tuy nhiên, vô cùng khó xảy ra trường hợp tất cả các sai số làm tròn chỉ về cùng một hướng. Nếu thay vào đó giả định rằng các sai số làm tròn là các biến làm tròn độc lập, thì tổng sai số làm tròn tỷ lệ thuận với $\backslash varepsilon\; /\; \backslash sqrt\{h\}$.

Vì vậy, đối với những giá trị kích thước bước cực nhỏ, sai số cắt cụt sẽ nhỏ nhưng tác động của sai số làm tròn có thể lớn. Hầu hết các tác động của sai số làm tròn có thể dễ dàng tránh được nếu phép tổng đền bù (compensated summation) được sử dụng trong việc xây dựng công thức cho phương pháp Euler.

Hiệu chỉnh và mở rộng

Một hiệu chỉnh đơn giản của phương pháp Euler loại bỏ các vấn đề ổn định đã lưu ý trong phần trước là phương pháp Euler lùi lại (backward): :$y\_\{n+1\}\; =\; y$n + h f(t{n+1}, y{n+1}). Phương pháp này khác với phương pháp Euler (tiêu chuẩn, hoặc tiếp tới) là hàm $f$ được đánh giá tại điểm cuối của bước, thay vì điểm xuất phát. Phương pháp Euler lùi lại là một phương pháp ẩn, có nghĩa là công thức của phương pháp Euler lùi lại có $y${n+1} ở cả hai bên, vì vậy khi áp dụng phương pháp Euler lùi lại chúng ta phải giải một phương trình. Điều này làm cho việc thực hiện tốn kém (thời gian,...) hơn.

Những hiệu chỉnh khác của phương pháp Euler đối với sự ổn định đã đưa đến phương pháp Euler mũ hoặc phương pháp Euler bán ẩn.

Các phương pháp phức tạp hơn có thể đạt được bậc cao hơn (và chính xác hơn). Một khả năng đó là sử dụng nhiều hơn các đánh giá hàm. Điều này được minh họa bằng phương pháp điểm giữa đã được đề cập trong bài viết này: :$y\_\{n+1\}\; =\; y\_n\; +\; h\; f\; \backslash Big(t\_n\; +\; \backslash tfrac12\; h,\; y\_n\; +\; \backslash tfrac12\; h\; f(t\_n,\; y\_n)\; \backslash Big).$ Điều này dẫn đến họ của các phương pháp Runge-Kutta.

Một khả năng khác là sử dụng nhiều hơn các giá trị quá khứ, như được minh họa bằng phương pháp Adams-Bashforth hai bước:  :$y\_\{n+1\}\; =\; y$n + \tfrac32 h f(t{n}, y{n}) - \tfrac12 h f(t{n-1}, y_{n-1}). Điều này dẫn họ của các phương pháp đa bước tuyến tính.