Chương 3: Tích phân của các PVT (Phần 1)

Xem tiếp Phần 2 * Trở về Mục lục cuốn sách

Trong chương này, ta sẽ xét những kĩ thuật số trị thông dụng phục vụ cho việc lấy tích phân của hệ các phương trình vi phân thường (PVT). Sau đó ta sẽ dùng những kĩ thuật này để mô phỏng quỹ đạo của những đường bay khác nhau của quả bóng chày.

Theo định nghĩa, một phương trình vi phân thường, hay PVT, là một phương trình vi phân trong đó tất cả các biến phụ thuộc đều là hàm của một biến độc lập đơn lẻ. Còn PVT bậc n là phương trình sao cho khi tối giản, đạo hàm bậc cao nhất mà nó chứa là bậc n.

Theo định luật Newton về chuyển động, chuyển động của một tập hợp vật rắn có thể được rút gọn về một hệ các PVT bậc 2 trong đó thời gian t là biến độc lập chung. Chẳng hạn, các phương trình chuyển động của hệ n chất điểm tương tác nhau trong không gian một chiều có dạng:
d2(xj) / dt2 = Fj(x1, . . . , xn, t) / mj
với j = 1 đến n, trong đó xj là vị trí của chất điểm thứ j, mj là khối lượng của nó, v.v. Lưu ý rằng một hệ n PVT bậc 2 luôn có thể viết lại thành một hệ gồm 2n PVT bậc nhất. Vì vậy, hệ phương trình chuyển động trên có thể được viết lại thành:
dxj / dt = vj, 
dvj / dt = Fj(x1, . . . , xn, t) / mj. 
với j = 1 đến n. Ta thấy được kiến thức giải bằng số trị hệ PVT bậc nhất có thể giúp cho việc quan sát động thái của nhiều hệ động lực thú vị.

Phương pháp Euler

Xét hệ PVT bậc nhất dưới dạng tổng quát sau (1.19)
yʹ = f(x, y), 
trong đó ʹ là để chỉ d / dx, chịu ảnh hưởng của điều kiện đầu dạng
y(x0) = y0. 
Rõ ràng là nếu ta có thể tìm được phương pháp số trị để giải bài toán này thì sẽ không khó để khái quát cách giải cho hệ n PVT bậc nhất.

Điều quan trọng là thấy được nghiệm số trị của một phương trình vi phân chỉ là giá trị xấp xỉ của nghiệm thực. Nghiệm thực, y(x), của PT (1.19) có thể là một hàm liên tục của biến x. Tuy nhiên, khi giải phương trình theo cách số trị thì điều tốt nhất chúng ta làm được chỉ là tính các giá trị xấp xỉ của hàm số y(x) tại một loạt các điểm nút rời rạc, xn trong đón = 0, 1, 2, ⋯x0 < x1 < x2. Ở đây, ta chỉ xét đến các điểm nút cách đều nhau, với
xn = x0 + n h. 
trong đó, đại lượng h được gọi là độ dài bước. Đặt yn là giá trị xấp xỉ của y(x) tại điểm nút xn. Một sơ đồ tích phân số trị chính là một phương pháp mà bằng cách nào đó dùng những thông tin có ở PVT gốc, PT (1.19), để lập ra một loạt quy tắc liên hệ các giá trị yn với nhau.

Sơ đồ tích phân đơn giản nhất được Leonhard Euler, nhà toán học Thụy Sĩ vĩ đại của thế kỉ thứ 18, tìm ra và do đó được gọi là phương pháp Euler. Cần nói thêm là hầu như tất cả những phương pháp tiêu chuẩn dùng trong phân tích số trị đều được tìm ra trước phát minh máy tính. Ngày xưa, nhà khoa học phải tính tay—hẳn là một công việc dài và tẻ nhạt! Giả sử rằng ta đã tính được một xấp xỉ, yn cho nghiệm y(x), của PT (e1.19) tại điểm nút xn. Từ đó, độ dốc xấp xỉ của y(x) tại điểm này được cho bởi
ynʹ = f(xn, yn). 
Ta hãy xấp xỉ đường cong y(x) như một đoạn thẳng nối hai điểm nút lân cận xnxn + 1. Điều này dẫn đến
yn + 1 = yn + ynʹ h, 
hay
yn + 1 = yn + f(xn, yn) h. 
Công thức trên là cơ sở cho phương pháp Euler. Nó cho phép ta tính tất cả giá trị yn khi biết trước giá trị đầu y0 tại điểm nút thứ nhất x0. Phương pháp Euler được minh họa ở Hình 1.

euler
Hình 1. Minh họa cho phương pháp Euler.

Sai số số trị

Có hai nguồn chính gây sai số trong sơ đồ tích phân số trị cho PVT, đó là sai số chặt cụtsai số làm tròn. Sai số chặt cụt nảy sinh trongphương pháp Euler vì đường cong y(x) nói chung không phải là một đoạn thẳng nối giữa hai điểm nút lân cận xnxn + 1, như đã giả sử. Sai số đi đôi với sự xấp xỉ này có thể dễ dàng tính được bằng cách khai triển Taylor cho hàm y(x) quanh điểm x = xn:
y(xn + h) = y(xn) + h yʹ(xn) + h2 / 2 yʹʹ(xn) + ⋯

= yn + h f(xn, yn) + h2 / 2 yʹʹ(xn) + ⋯.  (1.25)
So sánh các PT (1.24) và (1.25) ta có
yn + 1 = yn + h f(xn, yn) + O(h2). 
Nói cách khác, mỗi khi tính một bước bằng phương pháp Euler ta phạm phải sai số chặt cụt bằng O(h2), trong đó h là độ dài bước. Nếu giả sử ta dùng phương pháp Euler để lấy tích phân PVT trên một khoảng x với độ dài cỡ bằng 1 đơn vị.Điều này đòi hỏi phải có O(h - 1) bước. Nếu mỗi bước chứa một sai số O(h2), và các sai số cộng dồn với nhau (giả sử đơn giản này thiên về hướng “an toàn”), thì sai số tổng hợp sẽ là O(h). Nói cách khác, sai số đi kèm vớitích phân một PVT trên một khoảng hữu hạn bằng phương pháp Euler sẽ tỉ lệ thuận với độ dài bước. Vì vậy, nếu muốn giữ sai số tương đối của phép tích phân dưới mức 10 - 6, ta sẽ cần phải thực hiện khoảng 1 triệu điểm cho mỗi đơn vị độ dài lấy tích phân của x. Cần nói thêm là, phương pháp Euler thuộc loại phương pháp tích phân bậc nhất vì sai số chặt cụt gắn với việc tích phân trên khoảng hữu hạn có độ lớn cỡ h1. Nói rộng ra, một phương pháp tích phân được gọi là có bậc n nếu sai số chặt cụt của nó trong mỗi bướcO(hn + 1).

Lưu ý rằng sai số chặt cụt có thể xảy ra ngay cả khi máy tính thực hiện hoàn toàn chính xác phép tính với số có phần thập phân. Không may là máy tính không thể là như vậy. Thật ra, máy chỉ có thể lưu trữ được số có phần thập phân đến một lượng nhất định những chữ số sau dấu phẩy. Mỗi loại máy tính đều có một số đặc trưng, η, được định nghĩa là số nhỏ nhất sao cho khi cộng vào một số có độ lớn cỡ đơn vị thì sẽ cho ra một số mới. Mỗi phép tính với số có phần thập phân sẽ gây ra một sai số làm tròn bằng O(η) nảy sinh từ độ chính xác có hạn mà số có phần thập phân được lưu trong máy tính. Giả sử ta dùng phương pháp Euler để lấy tích phân PVT trên một khoảng x có độ dài cỡ đơn vị. Việc này cần đến O(h–1) bước tính tích phân, và do đó, O(h–1) phép tính với số có phần thập phân. Nếu mỗi phép tính như vậy gây ra lỗi bằng O(η), và các sai số có tính cộng dồn thì sai số làm tròn tổng hợp sẽ là O(η/h).

Sai số tổng hợp ε trong việc lấy tích phân PVT trên một khoảng x có độ dài cỡ đơn vị thì (xấp xỉ) bằng tổng của hai loại sai số chặt cụt và làm tròn. Do đó, với phương pháp Euler,
ε ~ η/h + h.
Rõ ràng là khi bước tính còn lớn, sai số chủ yếu sẽ là dạng chặt cụt, còn khi bước tính nhỏ thì sai số làm tròn là chủ yếu. Sai số tổng hợp nhận giá trị nhỏ nhất, ε0 ∼ η1/2, khi h = h0 ∼ η1/2. Rõ ràng là chẳng có ích gì khi co ngắn bước tính h xuống nhỏ hơn h0, vì điều này làm tăng số phép tính nhưng không làm tăng độ chính xác cuối cùng. Cũng thấy rõ rằng độ chính xác cuối cùng của phương pháp Euler (hay bất kì phương pháp tích phân nào khác) đều được quyết định bởi η, liên quan đến việc lưu trữ các số có phần thập phân trên máy.

Giá trị của η phụ thuộc vào số byte mà phần cứng máy tính dùng để lưu trữ số có phần thập phân. Với các máy tính dòng IBM-PC, giá trị thích hợp cho số thập phân có độ chính xác képη = 2,22 × 10–16 (giá trị này được nêu trong header file hệ thống, float.h). Điều đó dẫn đến là độ dài bước tối thiểu nên dùng cho phương pháp Euler khi chạy trên máy tính là h0 ∼ 10–8, ứng với sai số tích phân tương đối nhỏ nhất bằng ε0 ∼ 10–8. Mức chính xác này quá đủ đối với hầu hết các tính toán khoa học. Tuy nhiên, cần lưu ý rằng giá trị η tương ứng với số thập phân có độ chính xác đơn chỉ là η = 1,19 × 10–7, ứng với độ dài bước nên chọn bằng h0 ∼ 3 × 10–4 và sai số tương đối tối thiểu ε0 ∼ 3 × 10–4 khi dùng phương pháp Euler. Mức chính xác này nói chung không đáp ứng được yêu cầu tính toán khoa học, điều này lí giải việc các tính toán khoa học luôn được thực hiện theo số thực kép (double) thay vì đơn (single) trên các máy dòng IBM-PC (cũng như với hầu hết các dòng máy khác).

Bất ổn định số trị

Ta xét ví dụ sau đây. Giả sử rằng PVT là
yʹ =  - α y, 
trong đó α > 0, chịu ảnh hưởng của điều kiện biên
y(0) = 1. 
Dĩ nhiên là ta có thể giải bài toán này theo cách giải tích để tìm ra
y(x) = exp( - α x). 
Chú ý rằng nghiệm này là một hàm đơn điệu giảm theo x. Ta cũng có thể giải bài toán này theo cách số trị bằng phương pháp Euler. Các điểm nút thích hợp sẽ là
xn = n h, 
trong đó n = 0, 1, 2, ⋯. Phương pháp Euler cho ta
yn + 1 = (1 - α h) yn. 
Có một chi tiết đáng chú ý sau đây. Nếu h > 2 / α thì yn + 1∣ > ∣yn. Nói cách khác, nếu độ dài bước quá lớn thì nghiệm số trị sẽ là một hàm dao động theo x có biên độ đơn điệu tăng, tức là nghiệm số trị phân kì khỏi nghiệm thực. Dạng thất bại này của sơ đồ tích phân số trị được gọi là bất ổn định số trị. Tất cả sơ đồ tích phân đơn giản đều trở nên bất ổn định khi bước tính bị kéo dài quá mức.

Phương pháp Runge-Kutta

Có hai lý do chính khiến cho phương pháp Euler thường không được dùng trong tính toán khoa học. Thứ nhất, sai số do chặt cụt trong từng bước tính lớn hơn rất nhiều so với các phương pháp khác cao cấp hơn (với cùng giá trị của h). Thứ hai, phương pháp Euler rất dễ mắc phải sự bất ổn định số trị.

Những phương pháp thường được nhà khoa học dùng nhất để tính tích phân cho PVT được nghiên cứu phát triển đầu tiên bởi hai nhà toán học người Đức, C.D.T. Runge và M.W. Kutta vào nửa sau thế kỉ 19.1 Nguyên lý cơ bản của phương pháp Runge-Kutta sẽ được giới thiệu ngay sau đây.

Lí do chính khiến cho phương pháp Euler mắc sai số chặt cụt quá lớn trong mỗi bước tính là vì khi dẫn nghiệm từ xn đến xn + 1 phương pháp chỉ ước lượngđạo hàm tại điểm đầu đoạn, tức là tại xn. Vì vậy, phương pháp này rất không đối xứng giữa các điểm đầu và cuối đoạn. Ta có thể lập một phương pháp lấy tích phân đối xứng hơn bằng cách thêm vào một bước thử kiểu như Euler tại điểm giữa của đoạn, rồi dùng hai giá trị của cả xy tại điểm giữa này để làm thành bước thật cho cả đoạn. Cụ thể hơn,
k1 = h f(xn, yn), 

k2 = h f(xn + h/2 , yn + k1/2),

yn + 1 = yn + k2 + O(h3).
Như đã chỉ ra ở số hạng sai số, sự đối xứng này sẽ triệt tiêu sai số bậc nhất, làm cho phương pháp trở thành bậc hai. Thật ra, phương pháp trên thường được biết đến với tên gọi phương pháp Runge-Kutta bậc hai. Phương pháp Euler có thể được coi là phương pháp Runge-Kutta bậc nhất.

Dĩ nhiên, ta chẳng có lý do gì phải dừng lại ở phương pháp bậc hai. Bằng cách dùng hai bước tính thử trong mỗi khoảng, ta có thể triệt tiêu cả hai thành phần sai số bậc nhất và bậc hai, và do đó tạo nên một phương pháp Runge-Kutta bậc ba.Tương tự, với ba bước tính thử trong mỗi khoảng ta sẽ có phương pháp bậc bốn, va cứ như vậy.2

Biểu thức chung cho sai số tổng hợp, ε, phát sinh khi lấy tích phân PVT trên một khoảng x có độ dài cỡ đơn vị khi dùng phương pháp Runge-Kutta bậc n được tính xấp xỉ bằng
ε ~ η/h + hn.
Ở đây, số hạng thứ nhất tương ứng với sai số làm tròn, còn số hạng thứ hai biểu thị sai số chặt cụt. Bước tính nhỏ nhất nên chọn, h0, và sai số nhỏ nhất, ε0, nhận các giá trị
h0 ∼ η1 / (n + 1),
ε0 ∼ ηn / (n + 1). 

Ở Bảng 1, các giá trị này được liệt kê theo n với η = 2,22 × 10–16 (giá trị tương ứng với số có độ chính xác kép trên máy tính dòng IBM-PC). Có thể thấy rằng h0 tăng và ε0 giảm khi n tăng. Tuy nhiên, mức độ thay đổi tương đối của các đại lượng này ngày càng ít đi khi n tăng lên.

n h0 ε0
1 1. 5 × 10 - 8 1. 5 × 10 - 8
2 6. 1 × 10 - 6 3. 7 × 10 - 11
3 1. 2 × 10 - 4 1. 8 × 10 - 12
4 7. 4 × 10 - 4 3. 0 × 10 - 13
5 2. 4 × 10 - 3 9. 0 × 10 - 14

Bảng 1. Bước tính nhỏ nhất nên chọn, h0, và sai số nhỏ nhất, ε0, của phương pháp Runge-Kutta bậc n khi lấy tích phân trên một khoảng hữu hạn tính trên số có độ chính xác kép bằng máy tính dòng IBM-PC.

Trong phần lớn các trường hợp, yếu tố hạn chế khi ước tính tích phân số trị một PVT không phải là sai số làm tròn, mà là khối lượng tính toán cần thực hiện để tính hàm f(x, y). Lưu ý rằng, nhìn chung, một phương pháp Runge-Kutta bậc n yêu cần n lần ước lượng hàm này trong mỗi bước tính. Có thể thấy được là khi n tăng lên ta sẽ nhanh chóng đạt đến một giới hạn mà vượt khỏi nó, mọi ích lợithu được từ độ chính xác của phương pháp bậc cao hơn sẽ bị triệt tiêu bởi “phí tổn” do phải thực hiện thêm những lượng giá hàm f(x, y) ở mỗi bước. Mặc dù không có một quy tắc cố định và đơn giản để xác định giới hạn trên, song với hầu hết những bài toán trong vật lý tính toán giới hạn trên tương đương vớin = 4. Nói cách khác, trong nhiều tình huống phương pháp tích phân Runge-Kutta bậc bốn là một lựa chọn hợp lý xét đến sự cân bằng giữa yêu cầu giảm thiểu sai số chặt cụt và giảm thiểu khối lượng tính toán trong mỗi bước.

Phương pháp Runge-Kutta bậc bốn tiêu chuẩn có dạng sau: (1.38–1.42)
k1 = h f(xn, yn)
k2 = h f(xn + h/2, yn + k1/2)

k3 = h f(xn + h/2, yn + k2/2)
k4 = h f(xn + h, yn + k3)

yn+1 = yn + k1/6 + k2/3 + k3/3 + k4/6 + O(h5).
Chính phương pháp này sẽ được dùng trong suốt khóa học để lấy tích phân các PVT bậc nhất. Hi vọng rằng bạn sẽ dễ tìm thấy cách khái quát phương pháp này để giải hệ PVT bậc nhất.

Một ví dụ đoạn chương trình RK4 tính với bước cố định

Dưới đây là một đoạn chương trình ví dụ để tính tích phân Runge-Kutta bậc 4 (RK4) trong đó có dùng thư viện Blitz++ (xem Mục [mảng nhiều chiều, Ch2]).

// rk4_fixed.cpp  

/* 
  Function to advance set of coupled first-order o.d.e.s by single step 
  using fixed step-length fourth-order Runge-Kutta scheme 

     x        ... independent variable 
     y        ... array of dependent variables  
     h        ... fixed step-length 

  Requires right-hand side routine 

     void rhs_eval (double x, Array<double,1> y, Array<double,1>& dydx) 

  which evaluates derivatives of y (w.r.t. x) in array dydx 
*/ 

#include <blitz/array.h> 

using namespace blitz; 

void rk4_fixed (double& x, Array<double,1>& y,  
                void (*rhs_eval)(double, Array<double,1>, Array<double,1>&),  
                double h) 
{ 
  // Array y assumed to be of extent n, where n is no. of coupled 
  // equations 
  int n = y.extent(0); 

  // Declare local arrays 
  Array<double,1> k1(n), k2(n), k3(n), k4(n), f(n), dydx(n); 

  // Zeroth intermediate step  
  (*rhs_eval) (x, y, dydx); 
  for (int j = 0; j < n; j++) 
    { 
      k1(j) = h * dydx(j); 
      f(j) = y(j) + k1(j) / 2.; 
    } 

  // First intermediate step  
  (*rhs_eval) (x + h / 2., f, dydx); 
  for (int j = 0; j < n; j++) 
    { 
      k2(j) = h * dydx(j); 
      f(j) = y(j) + k2(j) / 2.; 
    } 

  // Second intermediate step 
  (*rhs_eval) (x + h / 2., f, dydx); 
  for (int j = 0; j < n; j++) 
    { 
      k3(j) = h * dydx(j); 
      f(j) = y(j) + k3(j); 
    } 

  // Third intermediate step  
  (*rhs_eval) (x + h, f, dydx); 
  for (int j = 0; j < n; j++) 
    { 
      k4(j) = h * dydx(j); 
    } 

  // Actual step  
  for (int j = 0; j < n; j++) 
    { 
      y(j) += k1(j) / 6. + k2(j) / 3. + k3(j) / 3. + k4(j) / 6.; 
    } 
  x += h; 

  return; 
}

Một ví dụ tính toán

Xét hệ PVT sau:
dx/dt = v
dv/dt = –k x,
với các điều kiện đầu x(0) = 0v(0) = √k tại t = 0. Thật ra hệ này có thể được giải theo cách giải tích để có nghiệm
x = sin(√k • t).
Hãy so sánh nghiệm trên với nghiệm được tính số trị bằng cả phương pháp Euler lẫn Runge-Kutta bậc 4. Hình 2 cho thấy các sai số trong quá trình tính tích phân theo hai phương pháp (với k = 1, từ t = 0 đến t = 10, rồi lấy hiệu số giữa các nghiệm số trị và giải tích) được vẽ lên biểu đồ loga theo bước tính h. Tất cả các phép tính đều được thực hiện với độ chính xác đơn, tức là dùng biến float thay vì double. Có thể thấy rằng với các giá trị h lớn, sai số của phương pháp Euler lớn hơn rất nhiều so với 1 (nghĩa là độ lớn của nghiệm số trị vượt quá nhiều so với nghiệm giải tích); điều này cho thấy sự bất ổn định số trị. Không có dấu hiệu bất ổn định tương tự trong phương pháp Runge-Kutta. Với h lớn vừa phải, sai số của phương pháp Euler giảm đều đặn với cỡ độ lớn h–1: trong vùng này, sai số chủ đạo là sai số chặt cụt, vốn có độ lớn cỡ h–1 với một phương pháp bậc nhất. Sai số của phương pháp Runge-Kutta có độ lớn cỡ h–4—như mong đợi với một sơ đồ bậc bốn— trong vùng mà sai số chặt cụt đóng vai trò chủ đạo. Lưu ý rằng khi h giảm, sai số gắn với cả hai phương pháp lại bắt đầu tăng theo một đường cong răng cưa với độ lớn xấp xỉ cỡ h1. Đây là biểu hiện của sai số làm tròn. Sai số nhỏ nhất của cả hai phương pháp ứng với biên giới giữa hai vùng lần lượt do sai số chặt cụt và sai số làm tròn chi phối. Như vậy với phương pháp Euler, sai số nhỏ nhất là khoảng 10–3 tại h ∼ 10–3, còn với phương pháp Runge-Kutta, sai số nhỏ nhất vào khoảng 10–5 tại h ∼ 10–1. Rõ ràng là phương pháp Runge-Kutta mạnh hơn nhiều so với phương pháp Euler, vì nó có khả năng đạt độ chính xác cao hơn rất nhiều mà chỉ dùng số bước tính ít hơn đáng kể (nghĩa là h lớn nhiều).

float
Hình 2. Sai số tổng hợp khi tính tích phân bằng phương pháp Euler (đường liền nét) và phương pháp Runge-Kutta bậc bốn (đường đứt nét) theo độ dài bước tính h. Tính toán thực hiện với độ chính xác đơn.

Hình 3 cũng cho thấy số liệu tương tự như trên Hình 2, nhưng bây giờ việc tính toán được thực hiện với độ chính xác kép. Hình vẽ thể hiện những nét khái quát như đã thấy ở Hình 2. Sự khác biệt chủ yếu là sai số làm tròn đã được giảm xuống thêm khoảng 9 cấp độ lớn; điều này cho phép phương pháp Runge-Kutta đạt được sai số nhỏ nhất vào tầm 10–12 (xem Bảng 1)—một sự thể hiện đáng nể!

double
Hình 3. Sai số tổng hợp khi tính tích phân bằng phương pháp Euler (đường liền nét) và phương pháp Runge-Kutta bậc bốn (đường đứt nét) theo độ dài bước tính h. Tính toán thực hiện với độ chính xác kép.

Hai hình 2 và 3 cho thấy tại sao nhà khoa học hiếm khi dùng phương pháp Euler, hoặc tính với số có độ chính xác đơn, để tìm tích phân của hệ PVT.

Các phương pháp tích phân thích ứng

Xét hệ PVT sau:
dx / dt = u,
du / dt = u / t - 4 k t2 x, 

với các điều kiện biên x = √k • ν2u = 2 √k • ν tại x = ν, trong đó 0 < ν ≪ 1. Hệ này có thể được giải theo cách giải tích để tìm ra
x = sin(√k • t2).
Một đặc điểm riêng của lời giải trên là cỡ độ dài biến đổi của nó giảm đi nhanh chóng khi t tăng.

fixed
Hình 4. Sai số tổng hợp của phép tính tích phân với bước tính cố định (h = 0,01), phương pháp Runge-Kutta bậc bốn, theo biến độc lập t, để giải hệ PVT trong đó cỡ độ dài biến đổi của nó giảm đi nhanh chóng khi t tăng. Tính toán thực hiện với độ chính xác kép.

Ta hãy so sánh nghiệm trên với nghiệm thu được bằng phương pháp Runge-Kutta bậc bốn. Hình 4 cho thấy sai số khi dùng phương pháp này (khi giải hệ trên với k = 10, từ t = 10–3 đến t = t, và sau đó lấy hiệu giữa nghiệm số trị và nghiệm giải tích) với bước tính cố định bằng h = 0,01, dưới dạng biểu đồ theo biến độc lập, t. Có thể thấy rằng, mặc dù ban đầu sai số rất nhỏ, nhưng nó tăng nhanh chóng khi độ dài cỡ biến đổi của nghiệm giảm đi (tức là khi t tăng), và nhanh chóng lớn đến mức không chấp nhận được. Dĩ nhiên, ta có thể giảm sai số bằng cách rút ngắn độ dài bước, h. Tuy nhiên, cách này rất không hiệu quả. Độ dài bước chỉ cần rút ngắn khi t lớn. Còn lúc t nhỏ thì không cần rút ngắn chút nào. Rõ ràng, giải pháp lý tưởng cho bài toán này là một phương pháp tích phân trong đó bước tính thay đổi được nhằm duy trì một sai số chặt cụt tương đối đồng đều qua mỗi bước tính. Một phương pháp tích phân thích ứng như vậy sẽ chọn bước tính dài khi sự thay đổi chiều dài kích cỡ còn lớn, và ngược lại.

Hãy tìm hiểu xem bằng cách nào ta có thể chuyển đổi phương pháp Runge-Kutta bậc bốn với chiều dài bước cố định sang phương pháp tương đương có tính thích ứng. Trước hết, ta cần ước tính được sai số chặt cụt ở mỗi bước. Giả sử rằng chiều dài bước hiện tại là h. Ta có thể ước tính sai số chặt cụt, ε của bước tính hiện thời bằng cách lấy hiệu giữa nghiệm thu được bằng bước qua h/2 hai lần và qua h một lần (bắt đầu từ cùng một điểm trong cả hai trường hợp). Gọi ε0 là sai số chặt cụt muốn có ở mỗi bước. Làm thế nào để điều chỉnh h sao cho đảm bảo được sai số chặt cụt trong bước tính tiếp theo không vượt quá giá trị trên? Ta thấy được, từ PT (1.42), sai số chặt cụt trong mỗi bước tính của sơ đồ bậc bốn tỉ lệ với h5. Vì vậy công thức điều chỉnh độ dài bước tính cần có dạng 3

h_{\rm new} = h_{\rm old}\left|\frac{\epsilon_0}{\epsilon}\right|^{1/5}.
Theo công thức này, độ dài bước tính phải tăng lên nếu sai số chặt cụt trong từng bước quá nhỏ, và ngược lại, sao cho sai số trong mỗi bước gần như giữ nguyên ở mức ε0.

Có một loạt những điều cần lưu ý trong nội dung trên. Với hệ gồm n PVT, sai số chặt cụt trong mỗi bước, ε, tất nhiên cần phải là một trung bình trọng số hợp lý của các sai số ở từng phương trình. Cũng có một câu hỏi đặt ra là liệu ε nên là một sai số tuyệt đối hay sai số tương đối. Sai số tương đối của phương trình thứ i đơn giản bằng sai số tuyệt đối chia cho yi, trong đó yi là giá trị hiện thời của biến phụ thuộc thứ i. Việc dùng sai số tuyệt đối thích hợp hơn trong trường hợp hệ phương trình mà biên độ của các biến phụ thuộc bị chặn. Việc dùng sai số tương đối thích hợp hơn trong trường hợp biên độ của biến phụ thuộc có chỗ bùng nổ, nhưng dấu của chúng vẫn giữ nguyên. Sau cùng, một ước lượng sai số hỗn hợp— thường là số nhỏ hơn giữa sai số và tuyệt đối—sẽ phù hợp với hệ trong đó biên độ của biến phụ thuộc có điểm bùng nổ và dấu của chúng dao động. Một cách làm hay là đặt ra giới hạn cho phép co giãn độ dài bước tính trong quá trình tính toán, chẳng hạn bằng cách không cho độ dài bước tăng hoặc giảm quá một hệ số S > 1 trong mỗi bước. Điều này ngăn cản không cho h dao động quá mức quanh giá trị tối ưu của nó. Hiển nhiên là nếu h quá nhỏ thì phương pháp tính tích phân sẽ thất bại, và cần phải dừng việc tính toán với một thông báo lỗi. Sau cùng, cần đặt ra một giới hạn quyết định xem h có thể lớn đến đâu—thật không may là các phương pháp thích ứng có xu hướng trở nên “lạc quan” thái quá khi việc tính tích phân được thuận lợi.

adaptive
Hình 5. Sai số tổng hợp trong cách tính tích phân bằng phương pháp Runge-Kutta bậc bốn [với bước tính cố định h = 0. 01] (nét liền) và phương pháp thích ứng tương đương (ε0 = 10–8)(nét chấm), được vẽ theo biến độc lập, t, cho hệ PVT trong đó sự thay đổi độ dài kích cỡ giảm đi nhanh chóng khi t tăng. Phép tính được thực hiện với độ chính xác kép.

Hình 5 cho thấy sai số tính tích phân của phương pháp Runge-Kutta bậc bốn với độ dài bước không đổi và phương pháp thích ứng tương đương—được thiết lập như đã xét ở trên—như hai hàm số của biến độc lập, t. Những sai số này được tính bằng cách lấy tích phân của hệ đang xét, với k = 10, từ t = 10–3 đến t = t, và sau đó tính hiệu giữa nghiệm số trị và nghiệm giải tích. Độ dài bước cố định trong phương pháp thứ nhất là h = 0,01. Sai số chặt cụt mong muốn với mõi bước tính trong phương pháp thứ hai được chọn là ε0 = 10–8. Có thể thấy rằng phương pháp thích ứng thể hiện tốt hơn nhiều so với phương pháp bước cố định, vì đã duy trì được sai số tích phân gần như đồng đều khi sự thay đổi bước kích cỡ của nghiệm giảm xuống, tức là t tăng lên. Hình {f6x} minh họa cách đạt được điều này. Hình vẽ cho thấy độ dài bước, h, của phương pháp thích ứng như một hàm của t. Nó còn cho thấy phương pháp thích ứng duy trì được sai số chặt cụt khá ổn định qua từng bước tính, bằng việc giảm h khi t tăng.

h
Hình 6. Độ dài bước tính h trong phương pháp tích phân thích ứng ở hình vẽ trước, ở đây được vẽ như hàm số theo biến độc lập t. Tính toán được thực hiện với độ chính xác kép.

Một đoạn chương trình ví dụ cho phương pháp RK4 bước thích ứng

Dưới đây là một ví dụ chương trình để tính phương pháp RK4 có bước thích ứng, dựa trên đoạn chương trình tính với bước cố định ở trước đó. Lưu ý rằng đoạn chương trình này tính lại tất cả những bước mà sai số chặt cụt vượt quá giá trị mong muốn acc.

// rk4_adaptive.cpp  

/* 
  Function to advance set of coupled first-order o.d.e.s by single step 
  using adaptive fourth-order Runge-Kutta scheme 

     x       ... independent variable 
     y       ... array of dependent variables 
     h       ... step-length 
     t_err   ... actual truncation error per step  
     acc     ... desired truncation error per step 
     S       ... step-length cannot change by more than this factor from 
                  step to step 
     rept    ... number of step recalculations           
     maxrept ... maximum allowable number of step recalculations           
     h_min   ... minimum allowable step-length 
     h_max   ... maximum allowable step-length 
     flag    ... controls manner in which truncation error is calculated     

  Requires right-hand side routine  

        void rhs_eval (double x, Array<double,1> y, Array<double,1>& dydx) 

  which evaluates derivatives of y (w.r.t. x) in array dydx. 

  Function advances equations by single step whilst attempting to maintain  
  constant truncation error per step of acc: 

    flag = 0 ... error is absolute 
    flag = 1 ... error is relative 
    flag = 2 ... error is mixed 

  If step-length falls below h_min then routine aborts 
*/ 

#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#include <blitz/array.h> 

using namespace blitz; 

void rk4_fixed (double&, Array<double,1>&,  
                void (*)(double, Array<double,1>, Array<double,1>&),  
                double); 

void rk4_adaptive (double& x, Array<double,1>& y,  
                   void (*rhs_eval)(double, Array<double,1>, Array<double,1>&), 
                   double& h, double& t_err, double acc,  
                   double S, int& rept, int maxrept,  
                   double h_min, double h_max, int flag) 
{ 
  // Array y assumed to be of extent n,  where n is no. of coupled 
  // equations 
  int n = y.extent(0); 

  // Declare local arrays 
  Array<double,1> y0(n), y1(n); 

  // Declare repetition counter 
  static int count = 0; 

  // Save initial data 
  double x0 = x; 
  y0 = y; 

  // Take full step  
  rk4_fixed (x, y, rhs_eval, h); 

  // Save data 
  y1 = y; 

  // Restore initial data  
  x = x0; 
  y = y0; 

  // Take two half-steps  
  rk4_fixed (x, y, rhs_eval, h/2.); 
  rk4_fixed (x, y, rhs_eval, h/2.); 

  // Calculate truncation error 
  t_err = 0.; 
  double err, err1, err2; 
  if (flag == 0) 
    { 
      // Use absolute truncation error  
      for (int i = 0; i < n; i++) 
        { 
          err = fabs (y(i) - y1(i)); 
          t_err = (err > t_err) ? err : t_err; 
        } 
    } 
  else if (flag == 1) 
    { 
      // Use relative truncation error 
      for (int i = 0; i < n; i++) 
        { 
          err = fabs ((y(i) - y1(i)) / y(i)); 
          t_err = (err > t_err) ? err : t_err; 
        } 
    } 
  else  
    { 
      // Use mixed truncation error  
      for (int i = 0; i < n; i++) 
        { 
          err1 = fabs ((y(i) - y1(i)) / y(i)); 
          err2 = fabs (y(i) - y1(i)); 
          err = (err1 < err2) ? err1 : err2; 
          t_err = (err > t_err) ? err : t_err; 
        } 
    } 

  // Prevent small truncation error from rounding to zero 
  if (t_err == 0.) t_err = 1.e-15; 

  // Calculate new step-length  
  double h_est = h * pow (fabs (acc / t_err), 0.2); 

  // Prevent step-length from changing by more than factor S 
  if (h_est / h > S) 
    h *= S; 
  else if (h_est / h < 1. / S) 
    h /= S; 
  else 
    h = h_est; 

  // Prevent step-length from exceeding h_max 
  h = (fabs(h) > h_max) ? h_max * h / fabs(h) : h; 

  // Abort if step-length falls below h_min 
  if (fabs(h) < h_min) 
    {  
      printf ("Error - |h| < hmin\n"); 
      exit (1); 
    } 

  // If truncation error acceptable take step  
  if ((t_err <= acc) || (rept >= maxrept)) 
    {   
      rept = count; 
      count = 0; 
    } 
  // If truncation error unacceptable repeat step  
  else  
    { 
      count++; 
      x = x0; 
      y = y0; 
      rk4_adaptive (x, y, rhs_eval, h, t_err, acc,  
                    S, rept, maxrept, h_min, h_max, flag); 
    } 

  return; 
}

Các phương pháp tính tích phân cao cấp

Dĩ nhiên, các phương pháp Runge-Kutta chưa phải là giải pháp cuối cùng để tính tích phân các PVT. Các phương pháp Runge-Kutta đôi khi còn được gọi là phương pháp một bước, vì chúng dẫn nghiệm từ xn đến xn + 1 mà không cần phải biết các nghiệm tại xn - 1, xn - 2, v.v. Còn một lớp rộng gồm các phương pháp tích phân phức tạp hơn, gọi là nhóm phương pháp nhiều bước, trong đó dùng các nghiệm đã tính ở xn - 1, xn - 2, v.v. để dẫn nghiệm từ xn đến xn + 1. Ví dụ cho loại này gồm có nhóm phương pháp Adams 4 và các phương pháp Ước tính–Hiệu chỉnh.5 Ưu điểm chính của các phương pháp Runge-Kutta là chúng dễ thực hiện, và rất ổn định, và có tính chất “tự khởi đầu”; nghĩa là, không như phương pháp nhiều bước, ta không cần phải xử lý một vài bước đầu tiên như trường hợp đặc biệt bằng cách tính tích phân một bước [để “mồi”]. Nhược điểm chính của các phương pháp Runge-Kutta là chúng đòi hỏi nhiều thời gian hơn đáng kể so với các phương pháp nhiều bước để đạt cùng độ chính xác, và chúng không dễ cho những ước lượng tổng thể tốt về sai số chặt cụt. Tuy nhiên, với những hệ động lực có cấu trúc đơn giản như sẽ xét trong khóa học này, lợi ích từ tính đơn giản và dễ dùng của các phương pháp Runge-Kutta trội hơn những nhược điểm về khối lượng tính toán của chúng.


  1. Numerical recipes in C: the art of scientific computing, W.H. Press, S.A. Teukolsky, W.T. Vettering, and B.R. Flannery (Cambridge University Press, Cambridge UK, 1992), p. 710.
  2. Handbook of mathematical functions, M. Abramowitzand I.A. Stegun (Dover, New York NY, 1965), p. 896.
  3. Numerical recipes in C: the art of scientific computing, W.H. Press, S.A. Teukolsky, W.T. Vettering, and B.R. Flannery (Cambridge University Press, Cambridge, England, 1992), p. 714.
  4. Numerical methods, R.W. Hornbeck (Prentice-Hall, Englewood Cliffs NJ, 1975), p. 196.
  5. ibid, p. 199.
Advertisements

2 phản hồi

Filed under Vật lý tính toán

2 responses to “Chương 3: Tích phân của các PVT (Phần 1)

  1. Pingback: Chương 1: Giới thiệu chung | Blog của Chiến

  2. Pingback: Phân Tích Pháp Số 1 Tối Nay | Vận May Bóng Đá

Trả lời

Mời bạn điền thông tin vào ô dưới đây hoặc kích vào một biểu tượng để đăng nhập:

WordPress.com Logo

Bạn đang bình luận bằng tài khoản WordPress.com Đăng xuất / Thay đổi )

Twitter picture

Bạn đang bình luận bằng tài khoản Twitter Đăng xuất / Thay đổi )

Facebook photo

Bạn đang bình luận bằng tài khoản Facebook Đăng xuất / Thay đổi )

Google+ photo

Bạn đang bình luận bằng tài khoản Google+ Đăng xuất / Thay đổi )

Connecting to %s