Chương 7: Phương trình sóng (Phần 1)

Trỏ về Mục lục cuốn sách

Giới thiệu chung

Phương trình sóng trong không gian một chiều có dạng

2ξ / ∂t2 = c2 (∂2ξ / ∂x2),    [wave1d]

Ta thường gặp trong vật lý nhiều đến nỗi không cần thiết phải liệt kê ví dụ về chúng. Ở đây, ξ thường là một dạng dịch chuyển hoặc nhiễu động, còn c vận tốc (không đổi) của sóng. Phương trình sóng có nghiệm hình thức sau

ξ(x, t) = F(x − ct) + G(x + ct), 

trong đó FG là các hàm tùy ý. Phương trình trên biểu thị những con sóng có hình dạng bất kì truyền với vận tốc c lần lượt theo các phương  +x −x, mà không thay đổi hình dạng.

Phương trình sóng, với dạng bậc hai theo cả không gian lẫn thời gian, có thể được viết thành hệ hai phương trình bậc nhất bằng cách định nghĩa các biến mới v = ∂ξ / ∂tθ =  − c ∂ξ / ∂x. Biểu diễn PT ([wave1d]) dưới dạng những biến mới này ta được

\begin{array}{ccc}  \partial v / \partial t + c\,\partial \theta / \partial x &=& 0,\\  \partial \theta / \partial t + c\, \partial v / \partial x &=& 0.\\  \end{array}

Lưu ý rằng khi giải phương trình sóng theo cách số trị, thường sẽ tiện hơn nếu ta viết nó dưới dạng một hệ hai phương trình bậc nhất như trên.

Phương trình chuyển tải 1 chiều

Phương trình sóng có liên hệ chặt chẽ đến phương trình chuyển tải, mà trong không gian 1 chiều nó có dạng

u / ∂t =  − v ∂u / ∂x.    [advection]

Phương trình này miêu tả sự chuyển tải thụ động của một trường vô hướng u(x, t) nào đó được dòng chảy đưa đi với vận tốc không đổi v. Vì phương trình chuyển tải có phần đơn giản hơn phương trình sóng nên ta sẽ xét đến nó trước. Phương trình chuyển tải có nghiệm hình thức sau

u(x, t) = F(x − vt), 

trong đó F là một hàm tùy ý. Nghiệm này mô tả một xung có hình dạng bất kì được dòng chảy có vận tốc không đổi v cuốn đi mà không thay đổi hình dạng.

Ta tìm nghiệm của PT ([advection]) trong miền xl ≤ x ≤ xh, tuân theo các điều kiện biên Dirichlet đơn giản u(xl) = u(xh) = 0. Như thường lệ, ta làm rời rạc theo thời gian trên lưới đều tn = t0 + nδt, với n = 0, 1, 2, ⋯. Hơn nữa, theo phương x, ta làm rời rạc trên lưới đều xi = xl + iδx, với i = 0, N + 1, trong đó δx = (xh − xl) / (N + 1). Khi áp dụng lược đồ sai phân hiện cho thời gian và một lược đồ sai phân trung tâm cho không gian, PT ([advection]) cho ta

(uin+1 − uin)  / δt =  − v (ui+1n − ui−1n)  / 2δx,     [ftcs1]

trong đó uin ≡ u(xi, tn). Phương trình trên có thể được viết lại thành

uin+1 = uin − (C / 2) (ui+1n − ui−1n),   [ftcs]

trong đó C = vδt / δx.

Ta hãy thực hiện một phân tích ổn định von Neumann cho lược đồ sai phân trên. Viết u(x, t) = û(t) exp(i kx), ta có ûin+1 = Aûin, trong đó

A = 1 − i C sin(kδx). 

Lưu ý rằng

A2 = 1 + C2 sin2(kδx) > 1. 

Như vậy, độ lớn của hệ số khuếch đại lớn hơn 1 với mọi k. Điều này có nghĩa là, thật không may, lược đồ sai phân đơn giản ([ftcs]) luôn bất ổn định không điều kiện.

Lược đồ Lax

Sự bất ổn định trong lược đồ sai phân ([ftcs]) có thể được khắc phục bằng cách thay thế uin ở vế phải bằng trung bình trong không gian của u lấy với các điểm nút lưới kề bên. Từ đó ta nhận được

uin+1 = ½ (ui+1n + ui−1n) − (C/2) (ui+1n − ui−1n),    [lax]

vốn được gọi là lược đồ Lax. Thực hiện phân tích ổn định von Neumann cho lược đồ Lax cho ra biểu thức độ khuếch đại sau:

A = cos(kδx) − i C sin(kδx). 

Đến đây

A2 = 1 − (1 − C2) sin2(kδx).     [laxa]

Hệ quả là lược đồ Lax ổn định không điều kiện, tức là A∣ < 1 với mọi k, chỉ cần C < 1. Từ định nghĩa của C, bất đẳng thức C < 1 có thể được viết thành

δt < δx / v. 

Đây là tiêu chuẩn ổn định Courant-Friedrichs-Lewy (hay CFL) nổi tiếng. Thực ra, tất cả lược đồ sai phân hiện ổn định để giải phương trình chuyển tải thì đều tuân theo điều kiện ràng buộc CFL, vốn quyết định đến bước thời gian dài nhất cho phép.

Dưới đây là đoạn chương trình giải phương trình chuyển tải một chiều bằng phương pháp Lax.

// Lax1D.cpp 

// Function to evolve advection equation in 1-d: 

//  du / dt + v du / dx = 0  for  xl <= x <= xh 

//  u = 0  at x=xl and x=xh 

// Array u assumed to be of extent N+2. 

// Now, ith element of array corresponds to 

//  x_i = xl + i * dx    i=0,N+1 

// Here, dx = (xh - xl) / (N+1) is grid spacing. 

// Function evolves u by single time-step. 

// C = v dt / dx, where dt is time-step. 

// Uses Lax scheme. 

#include <blitz/array.h> 

using namespace blitz; 

void Lax1D (Array<double,1>& u, double C) 
{ 
  // Set N. Declare local array. 
  int N = u.extent(0) - 2; 
  Array<double,1> u0(N+2); 

  // Evolve u 
  u0 = u; 
  for (int i = 1; i <= N; i++) 
    u(i) = 0.5 * (u0(i+1) + u0(i-1)) - 0.5 * C * (u0(i+1) - u0(i-1)); 

  // Set boundary conditions 
  u(0) = 0.; 
  u(N+1) = 0.; 
}

Hình [decay] cho thấy một ví dụ tính toán trong đó dùng chương trình trên để tính chuyển tải một xung Gauss. Điều kiện ban đầu là

u(x, 0) = exp[ − 100 (x − 0,5)2],      [initcn]

và tính toán được thực hiện với v = 1, δt = 2,49 × 10−3, và N = 200. Ngoài ra, xl = vtxh = 1 + vt. Lưu ý rằng C = 0, 5 với những tham số này. Có thể thấy được rằng xung này được đẩy đi với vận tốc đúng, nghĩa là xung xuất hiện gần như là tĩnh tại khi được vẽ theo trục x − vt. Không may là, xung không còn bảo toàn được hình dạng như lẽ ra phải vậy. Thay vào đó, xung trở nên dần thấp và bẹt xuống trong quá trình lan rồi, và cuối cùng tan biến hoàn toàn.

decay

Hình. Sự truyền tải của xung Gauss 1 chiều. Việc tính toán số trị được thực hiện với v = 1; δt = 2, 49 × 10−3, và N = 200. Đường nét liền biểu thị điều kiện ban đầu lúc t = 0, đường gạch ngắn là nghiệm số trị lúc t = 1, đường gạch dài là nghiệm số trị lúc t = 2, còn đường gạch chấm là nghiệm số trị lúc t = 3. [decay]

Rõ ràng là từ tính toán trên, có thể thấy lược đồ Lax giới thiệu một hiệu ứng giả vào trong bài toán chuyển tải. Ta có thể hiểu nguồn gốc của hiệu ứng này bằng cách tìm một nghiệm Fourier u(x, t) = ûk(t) exp( i kx), của PT ([advection]). Có thể dễ dàng tìm được

ûk(t) = ûk(0) exp( − i kvt). 

Lưu ý rằng ûk là một hằng số theo thời gian với mọi giá trị của k. Nói cách khác, biên độ của các hàm điều hòa Fourier trong nghiệm đúng của một phương trình chuyển tải vẫn không đổi theo thời gian—đó là pha của các hàm điều hòa trong quá trình phát triển. Bây giờ hãy xét PT ([laxa]). Có thể thấy rằng, miễn là điều kiện CFL C < 1 được thỏa mãn, độ lớn của hệ số khuếch đại, A, sẽ nhỏ hơn 1 với tất cả hàm điều hòa Fourier. Nói cách khác, lược đồ sai phân Lax khiến cho các hàm điều hòa Fourier tắt dần theo thời gian. Chính là sự tắt dần không đúng bản chất vật lý này của các hàm điều hòa Fourier đã gây ra hiệu ứng phân tán mạnh như được minh họa trên Hình [decay].

inst

Hình. Sự chuyển tải xung Gauss theo 1 chiều. Tính toán số trị được thực hiện với v = 1; δt = 9,95 × 1o−3, và N = 200. Đường chấm chấm biểu thị điều kiện ban đầu lúc t = 0,00; đường gạch ngắn là nghiệm số trị lúc t = 0,15; đường gạch dài là nghiệm số trị lúc t = 0,30; còn đường liền nét là nghiệm số trị lúc t = 0,37. [inst]

Hình [inst] cho thấy kết quả tính toán bằng lược đồ Lax trong trường hợp điều kiện CFL bị vi phạm. Tính toán này giống như lần trước, chỉ khác là bước thời gian được kéo dài tới δt = 9,95 × 10−3, cho ra số CFL, C = 2,0, vượt quá 1. Có thể thấy rằng biên độ xung tăng lên, và cuối cùng xung bị vỡ do mất ổn định khi bước sóng quá ngắn.

Lược đồ Crank-Nicolson

Lược đồ ẩn Crank-Nicolson để giải phương trình khuếch tán (xem Mục [Crank-Nicolson, Ch7]) có thể được điều chỉnh để giải phương trình khuếch tán. Theo đó, lấy trung bình vế phải của PT ([advection]) giữa các thời điểm đầu và cuối bước thời gian, ta thu được lược đồ sai phân dưới đây:

uin+1 + (C/4) (ui+1n+1 − ui−1n+1) = uin − (C/4) (ui+1n − ui−1n). 

Việc phân tích von Neumann của biểu đồ trên cho ta biểu thức sau về hệ số khuếch đại:

\displaystyle{  A = \frac{1 - \textrm{i} \,(C/2)\,\sin (k \,\delta x)}{1 + \textrm{i} \,(C/2)\,\sin (k \,\delta x)}.  }

Lưu ý rằng A∣ = 1 với mọi giá trị của k, bất kể giá trị của C. Điều đó nghĩa là lược đồ ẩn Crank-Nicolson không bị ràng buộc về số CFL, C < 1 (vì không có giá trị của k nào khiến A∣ > 1). Hơn nữa, không có sự tắt dần những thành phần điều hòa Fourier của nghiệm (vì A∣ = 1). Như vậy, không giống lược đồ Lax, ở lược đồ Crank-Nicolson ta không thấy xuất hiện sự phân tán mạnh do nguyên nhân số trị trong bài toán chuyển tải.

Dưới đây là một đoạn chương trình để giải phương trình chuyển tải 1 chiều bằng phương pháp Crank-Nicolson.

// Advect1D.cpp 

// Function to evolve advection equation in 1-d: 

//  du / dt + v du / dx = 0  for  xl <= x <= xh 

//  u = 0  at x=xl and x=xh 

// Array u assumed to be of extent N+2. 

// Now, ith element of array corresponds to 

//  x_i = xl + i * dx    i=0,N+1 

// Here, dx = (xh - xl) / (N+1) is grid spacing. 

// Function evolves u by single time-step. 

// C = v dt / dx, where dt is time-step. 

// Uses Crank-Nicolson scheme. 

#include <blitz/array.h> 

using namespace blitz; 

void Tridiagonal (Array<double,1> a, Array<double,1> b, Array<double,1> c,  
                  Array<double,1> w, Array<double,1>& u); 

void Advect1D (Array<double,1>& u, double C) 
{   
  // Find N. Declare local arrays. 
  int N = u.extent(0) - 2; 
  Array<double,1> a(N+2), b(N+2), c(N+2), w(N+2); 

  // Initialize tridiagonal matrix 
  for (int i = 2; i <= N;   i++) a(i) = - 0.25 * C; 
  for (int i = 1; i <= N; i++) b(i) = 1.; 
  for (int i = 1; i <= N-1; i++) c(i) = + 0.25 * C; 

  // Initialize right-hand side vector 
  for (int i = 1; i <= N; i++) 
    w(i) = u(i) - 0.25 * C * (u(i+1) - u(i-1)); 

  // Invert tridiagonal matrix equation 
  Tridiagonal (a, b, c, w, u); 

  // Calculate i=0 and i=N+1 values 
  u(0) = 0.; 
  u(N+1) = 0.; 
}

Hình [decay] cho thấy ví dụ tính bằng chương trình trên để giải bài toán chuyển tải một xung Gauss. Điều kiện ban đầu được chỉ định trong PT ([initcn]); tính toán được thực hiện với v = 1; δt = 9,95 × 10−3, và N = 200. Lưu ý rằng, với các tham số trên, C = 2,0. Có thể thấy rằng xung lan truyền với tốc độ (gần như là) đúng, và duy trì được hình dạng xấp xỉ với ban đầu. Rõ ràng là hiệu năng của lược đồ Crank-Nicolson mạnh hơn nhiều so với lược đồ Lax.

Hình. Sự chuyển tải của xung Gauss một chiều. Việc tính toán số trị được thực hiện với v = 1; δt = 9,95 × 10−3, và N = 200. Đường liền nét biểu thị điều kiện ban đầu lúc t = 0, đường gạch gắn là nghiệm số trị lúc t = 1, đường gạch dài là nghiệm số trị lúc t = 2, còn đường chấm gạch là nghiệm số trị lúc t = 3. [nodisperse]

Sai phân ngược dòng

Có thể sẽ thông cảm được, nếu ta sớm kết luận rằng lược đồ Crank-Nicolson là một phương pháp số đa năng, hiệu quả, và chính xác để giải phương trình chuyển tải. Phải công nhận rằng điều này đúng, nhưng chỉ trong trường hợp các hình sóng tương đối trơn tru. Thật không may là lược đồ Crank-Nicolson thực hiện rất kém khi tính chuyển tải dạng xung có nhánh lên hoặc nhánh xuống rất dốc. Điều này được thể hiện trên Hình [upwind1], trong đó cho thấy tính toán mà lược đồ Crank-Nicolson được dùng để chuyển tải một sóng hình vuông. Có thể thấy rằng những dao động không có thật được phát sinh ở cả đường lên và đường xuống của hình sóng. Hó ra rằng tất cả những lược đồ sai phân trung tâm để giải phương trình chuyển tải đều có vấn đề tương tự.

upwind1

Hình. Sự chuyển tải của một sóng vuông trong không gian 1 chiều. Tính toán số trị được thực hiện với v = 1, δt = 2,49 × 10−3, và N = 200. Đường chấm chấm biểu thị điều kiện ban đầu lúc t = 0,0, còn đường liền nét biểu thị nghiệm số trị lúc t = 0,1.[upwind1]

Cách duy nhất đã biết để loại trừ những dao động phi lý ở hai cạnh bên của một hình sóng biến đổi đột ngột là dùng một lược đồ sai phân ngược dòng. Ở lược đồ như vậy, các sai phan trong không gian bị kéo lệch đi về phía “ngược dòng”, nghĩa là phía mà dòng chảy khởi nguồn. Do vậy, dạng ngược dòng của lược đồ sai phân hiện đơn giản ([ftcs1]) là

(uin+1 − uin) / δt =  − v (uin − ui−1n) / δx, 

hay

uin+1 = uin − C (uin − ui−1n), 

Lưu ý rằng lược đồ này chỉ là bậc nhất theo không gian, trong khi tất cả các lược đồ khác đã đề cập đến đều là bậc hai. Phép phân tích ổn định von Neumann cho lược đồ trên có kết quả

A = 1 − C [1 − cos(kδx)] − i C sin(kδx). 

Lưu ý rằng

A2 = 1 − 2 C (1 − C) [1 − cos(kδx)]. 

Điều đó cho thấy A∣ < 1 với mọi k miễn là C < 1. Vì vậy, lược đồ sai phân ngược dòng là ổn định, chỉ cần điều kiện CFL được thỏa mãn. Hình [upwind2] biểu diễn tính toán trong đó lược đồ trên được dùng để mô phỏng sự chuyển tải một sóng hình vuông. Bây giờ không còn dao động phi lý ở hai bên của hình sóng. Mặt khác, hình sóng còn cho thấy biểu hiện của sự phân tán. Lược đồ sai phân ngược dòng thật sự cũng bị tình trạng phân tán phi lý giống như ở lược đồ Lax. Không may là hiện chưa có lược đồ sai phân đã biết nào vừa không bị phân tán vừa xử lý tốt trường hợp sóng có cạnh bên dựng đứng. Trên thực tế, những chương trình phức tạp để giải phương trình chuyển tải (hay sóng) nói chung đều dùng lược đồ ngược dòng ở những vùng sát với các cạnh bên dựng đứng của sóng (còn gọi là shock), và ở những chỗ khác thì dùng một lược đồ không phân tán có độ chính xác cao hơn.

upwind2

Hình. Sự chuyển tải của một sóng vuông trong không gian 1 chiều. Tính toán số trị được thực hiện với v = 1, δt = 2,49 × 10−3, và N = 200. Đường chấm chấm biểu thị điều kiện ban đầu lúc t = 0,0, còn đường liền nét biểu thị nghiệm số trị lúc t = 0,1.[upwind2]

Cần nói thêm là, ta dễ dàng thấy được rằng lược đồ sai phân xuôi dòng,

(uin+1 − uin) / δt =  − v (ui+1n − uin) / δx, 

không ổn định vô điều kiện.

Advertisements

1 phản hồi

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

One response to “Chương 7: Phương trình sóng (Phần 1)

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

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 Log Out / Thay đổi )

Twitter picture

Bạn đang bình luận bằng tài khoản Twitter Log Out / Thay đổi )

Facebook photo

Bạn đang bình luận bằng tài khoản Facebook Log Out / Thay đổi )

Google+ photo

Bạn đang bình luận bằng tài khoản Google+ Log Out / Thay đổi )

Connecting to %s