ODE and Runge-Kutta method

Bài này sẽ nói về maths. Tình hình là sắp tới sẽ còn phải cày bừa cái của nợ này nhiều.

1. Phương trình vi phân

Để nói tới ODE (Ordinary differential equation) thì trước tiên phải nói về phương trình vi phân (PTVP).

PTVP là một phương trình toán học nhằm biễu diễn mối quan hệ giữa một hàm chưa được biết (một hoặc nhiều biến) với đạo hàm của nó (có bậc khác nhau). Hàm chưa biết trong PTVP gọi là ẩn hàm. Mục tiêu của PTVP là tìm ra công thức của ẩn hàm, nhưng thông thường kết quả là một họ các hàm sai khác nhau một hằng số C. Hằng số này sẽ được xác định nếu có điều kiện ban đầu (initial value problem – IVP) hoặc điều kiện biên (boundary value problem – BVP) kèm theo.

Có các loại PTVP sau:

  • PTVP thường (Ordinary differential equationODE): là PTVP mà ẩn hàm là hàm một biến (độc lập). Lưu ý là ẩn hàm cũng có thể là hàm cho kết quả là vector, ma trận… chứ không nhất thiết phải là hàm có giá trị thực hay phức.
  • PTVP riêng phần (Partial differential equation – PDE): là PTVP mà ẩn hàm là hàm nhiều biến độc lập, và trong PTVP chứa các đạo hàm riêng của nó.
  • Delay differential equationDDE: PTVP mà ẩn hàm là hàm một biến, thường là thời gian. Đạo hàm của các ẩn hàm tại thời điểm nào đó được xác định dựa trên các thời điểm trước đó.
  • Stochastic differential equationSDE
  • Differential algebraic equationDAE

Mỗi loại này còn chia ra thành 2 loại là tuyến tính (linear) và phi tuyến (nonlinear). PTVP tuyến tính là các PTVP mà ẩn hàm và các đạo hàm của nó đều có lũy thừa 1, nếu có lũy thừa bậc cao hơn thì gọi là PTVP phi tuyến. Ví dụ với u là ẩn hàm thì u'=\alpha u là PTVP tuyến tính, nhưng u'=\alpha u^2 là PTVP phi tuyến.

Trong các PTVP tuyến tính, người ta còn quan tâm đặc biệt đến các PTVP tuyến tính có hệ số hằng (hệ số của ẩn hàm và các đạo hàm của nó trong PTVP đều là hằng). Đây là loại PTVP phức tạp nhất mà SV Đại học không chuyên ngành Toán ở VN được học cách giải.

2. Phương trình vi phân thường (ODE) và cách giải

2.1. Định nghĩa ODE

Gọi y là hàm theo biến x

y : R \rightarrow R

y^{(n)} là đạo hàm bậc n của y. Gọi F là hàm

F: R^{n+1} \rightarrow R

Khi đó phương trình có dạng

F(x, y, y', ..., y^{(n-1)}) = y^{(n)}

gọi là PTVP thường (ODE) bậc n.

PTVP có dạng như trên còn gọi là ODE rõ (explicit ODE) bậc n. Tổng quát hơn, PTVP có dạng

F(x, y, y', ..., y^{(n-1)}, y^{(n)}) = 0

được gọi là ODE ẩn (implicit ODE) bậc n, với F là hàm F: R^{n+2} \rightarrow R.

Theo định nghĩa trên, một PTVP gọi là tuyến tính (linear) nếu F có thể viết dưới dạng tổ hợp tuyến tính của các đạo hàm của ẩn hàm:

\displaystyle F(x, y, y', ... , y^{(n-1)}) = \sum_{i=0}^{n-1}{a_i(x)y^{(i)}} + r(x) = y^{(n)}

với a_i(x), r(x) là các hàm liên tục theo x. Nếu r(x) = 0 thì PTVP tuyến tính gọi là thuần nhất (homogeneous), ngược lại gọi là không thuần nhất (non-homogeneous, inhomogeneous).

Một ví dụ khá phổ biến về ODE có thể kể tới là phương trình định luật 2 Newton: F = ma, với a = \frac{d^2x}{dt^2}, ẩn hàm là x(t): [0, +\infty) \rightarrow R.

2.2. Cách giải

Một số ODE dạng đặc biệt có thể giải được với kết quả chính xác, tuy nhiên không phải ODE nào cũng giải được như vậy. Phần này trình bày cách giải một số ODE đặc biệt, phần cuối cùng sẽ trình bày thuật toán Runge-Kutta như một cách tìm lời giải xấp xỉ trong trường hợp không thể tìm được lời giải tổng quát.

2.2.1. PTVP thường tuyến tính bậc một phân tích được (Separable first-order linear ODE)

Dạng chung của PTVP loại này là

\displaystyle \frac{dy}{dx} + f(x)y = 0

với f(t) là hàm biết trước. Ta viết PT này thành

\displaystyle \frac{dy}{y} = -f(x)dx

Lấy nguyên hàm 2 vế được:

\displaystyle ln|y| = (\int {-f(x)\, dx}) + C

\Leftrightarrow y = \pm e^{(\int {-f(x)\, dx} + C)}

\Leftrightarrow y = \pm e^{ \int {-f(x)\, dx}}. e^C

\Leftrightarrow y = Ae^{ \int {-f(x)\, dx}}

Có thể kiểm chứng dễ dàng bằng cách thay y vào PTVP ban đầu.

2.2.2. PTVP thường tuyến tính bậc 1 không phân tích được (Non-separable first-order linear ODE)

Cho PTVP có dạng sau:

\displaystyle \frac{dy}{dx} + p(x)y = q(x)

Để giải dạng này, ta chọn một hàm làm thừa số tích phân:

\mu = e^{\int {p(x)dx}}

\mu có đặc điểm là đạo hàm của nó bằng chính nó nhân với p(x):

\displaystyle \frac{d\mu}{dx}= e^{\int {p(x)dx}}. p(x) = \mu p(x)

Nhân cả 2 vế của PTVP với \mu:

\displaystyle \mu  \frac{dy}{dx} + \mu p(x)y = \mu q(x)

Do \mu p(x) = \frac{d\mu}{dx} nên pt tương đương:

\displaystyle \mu  \frac{dy}{dx} + y\frac{d\mu}{dx} = \mu q(x)

\displaystyle \Leftrightarrow \frac{d}{dx}(\mu y) = \mu q(x)

Nguyên hàm 2 vế:

\displaystyle \mu y = (\int{\mu q(x) dx}) + C

\displaystyle \Leftrightarrow y = \frac{\int{\mu q(x) dx} + C}{\mu}

Do \mu cũng là hàm theo x nên vế phải có thể rút gọn được.

2.2.3. Cách giải PTVP thường tổng quát

Rất nhiều PTVP không rời vào các dạng đặc biệt để có thể tìm được lời giải một cách chính xác. Trong những trường hợp đó, người ta thường phải chấp nhận lời giải xấp xỉ bằng các phương pháp số học (gọi là Numerical ODEs), hoặc sử dụng giải tích để viết các chuỗi khai triển của lời giải. Một số phương pháp số học hay dùng để giải xấp xỉ PTVP bậc 1 là phương pháp Euler, Euler ngược… các phương pháp này là tiền đề của các phương pháp khác cho ra kết quả chính xác hơn như Runge-Kutta, Richardson extrapolation

Phần kế tiếp sẽ nói về phương pháp Runge-Kutta.

3. Phương pháp Runge-Kutta bậc 4 (fourth-order Rungke-Kutta method)

Đây là phương pháp thông dụng nhất trong họ các phương pháp Runge-Kutta, thường được gọi là RK4, phương pháp Runge-Kutta cổ điển hay đơn giản là phương pháp Runge-Kutta.

Cho PTVP với điều kiện ban đầu như sau:

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

Theo phương pháp RK4 thì lời giải cho PT này là:

y_{n+1} = y_n + \frac{1}{6}h(k_1 + 2k_2 + 2k_3 + k_4)

t_{n+1} = t_n + h

trong đó y_{n+1} là xấp xỉ của y(t_{n+1}), h là khoảng cách giữa 2 lần xấp xỉ liên tiếp, và k_1, k_2, k_3, k_4 là các slope, lần lượt là các hệ số góc trong khoảng h:

k_1 = f(t_n, y_n)

k_2 = f(t_n + \frac{1}{2}h, y_n + \frac{1}{2}hk_1)

k_3 =  f(t_n + \frac{1}{2}h, y_n + \frac{1}{2}hk_2)

k_4 =   f(t_n + h, y_n + hk_3)

Tổng quát hơn, ta có phương pháp Runge-Kutta bậc s, khi đó các công thức tương ứng được viết:

\displaystyle y_{n+1} = y_n + h\sum_{i=1}^s{b_i k_i},

trong đó:

k_1 = f(t_n, y_n),

k_2 = f(t_n + c_2 h, y_n + a_{21}hk_1),

k_3 = f(t_n + c_3 h, y_n + a_{31}hk_1 + a_{32}hk_2),

k_s = f(t_n + c_s h, y_n + \sum_{j=1}^{s-1}{a_{sj}hk_j})

Chi tiết có thể tham khảo trên wiki.

4. Ví dụ

Xét chuyển động của chất điểm khối lượng m, kích thước không đáng kể. Theo định luật 2 Newton, ta có:

F = ma

trong đó F là hợp lực tác dụng lên vật, a = \frac{d^2x}{dt^2} là gia tốc, x là tọa độ của vật trên quỹ đạo. Trong trường hợp này, x là ẩn hàm cần tìm. Ta có:

F = m\frac{d^2x}{dt^2}

\Leftrightarrow \frac{d^2x}{dt^2} = \frac{F}{m}

Phương trình này xuất hiện đạo hàm bậc 2 của x. Ta đặt ẩn phụ để đưa PT này thành 2 PT chỉ chứa đạo hàm bậc 1:

\frac{dx}{dt} = z, (1)

\frac{dz}{dt} = \frac{F}{m} (2)

Ta sẽ tìm x theo 2 cách, một là dùng phương pháp chỉ ra trong phần 2.2.1 để tìm công thức chính xác của x, 2 là dùng RK4 để tìm các giá trị của x trong đoạn cho trước. Kết quả sẽ cho thấy RK4 cho ra lời giải xấp xỉ khá tốt so với công thức chính xác.

Để so sánh kết quả giữa 2 phương pháp, ta cần một số điều kiện biên: \frac{F}{m} = 4, x(0) = 1, x(1) = 4.

4.1. Tìm nghiệm chính xác

Từ (2) ta được:

z = \frac{F}{m}t + C_1

thay vào (1):

\frac{dx}{dt} = \frac{F}{m}t + C_1

\Rightarrow dx = \frac{F}{m}t dt + C_{1}dt

\Rightarrow x = \frac{1}{2}\frac{F}{m}t^2 + C_{1}t + C_2

Do \frac{F}{m} = 4, x(0) = 1, x(1) = 4, ta có phương trình x theo t là:

x = 2t^2 + t + 1

Trong Hình 1 là đồ thị của x theo t trên đoạn [1, 3].

 

Hình 1: Đồ thị x(t) trên đoạn 1;3

4.2. Sử dụng RK4

Đoạn chương trình đơn giản sau cài đặt phương pháp RK4 để tìm x:

function w = runge_kutta4(a, b, N, alpha)

 h  = (b-a)/N;        %the step size
 t = zeros(1, N+1);
 w = zeros(1, N+1);
 t(1) = a;
 w(1) = alpha;        %the initial value

 for i = 1:N
  k1 = h*f(t(i), w(i));
  k2 = h*f(t(i)+h/2, w(i)+(k1)/2);
  k3 = h*f(t(i)+h/2, w(i)+(k2)/2);
  k4 = h*f(t(i)+h, w(i)+k3);
  w(i+1) = w(i) + (k1 + 2*k2 + 2*k3 + k4)/6;
  t(i+1) = a + i*h;
 end

 [t' w']

 plot(t, w, 'r*')

function dx = f(t, x)
 dx = 4*t + 1;

Ở đây hàm f được dùng để tính đạo hàm của dãy. Do x = 2t^2 + t + 1 nên dx=4t+1. Trong Matlab, thực thi lệnh sau:

>> runge_kutta4(1, 3, 2000, 4)

sẽ cho kết quả như Hình 2.

Hình 2: Đồ thị x(t) trên 1;3 tính theo phương pháp RK4

Trong matlab, phương pháp Runge-Kutta 4/5 được cài đặt sẵn trong hàm ode45, cùng với khá nhiều tùy chọn như tự động chọn step size, hỗ trợ event detection v.v… Ngoài ra matlab cũng hỗ trợ nhiều phương pháp khác để giải ODE.

của nó (có bậc khác nhau)
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s