Chương 5: Phương trình Poisson (Phần 1)

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

Giới thiệu chung

Trong chương này, ta sẽ tìm hiểu một số kỹ thuật số trị đơn giản để giải phương trình Poisson:

2u(r) = v(r). 

Ở đây, u(r) thường là một dạng trường thế, còn là một v(r) nguồn. Nghiệm của phương trình trên nói chung được tìm trong một thể tích V liên tục1 nào đó được bao kín bởi một mặt S. Có hai loại điều kiện biên chính cho phương trình Poisson. Trong điều kiện biên Dirichlet, thế u được chỉ ra trên mặt S. Trong điều kiện biên Neumann, gra-đien vuông góc với bề mặt của thế, u ⋅ dS được chỉ ra trên mặt bao.

Phương trình Poisson đặc biệt quan trọng trong tĩnh điện học và thuyết trọng trường Newton. Trong tĩnh điện học, ta có thể viết trường điện tích E dưới dạng một điện thế φ:

E =  − ∇ φ. 

Bản thân thế này đã thỏa mãn phương trình Poisson:

2φ =  − ρ / ε0,     [estat]

trong đó ρ(r) là mật độ điện tích, còn ε0 là độ điện thẩm của môi trường. Trong trường trọng lực Newton, ta có thể viết lực f tác dụng lên một đơn vị khối lượng theo thế trọng trường φ:

f =  − ∇ φ. 

Thế này thỏa mãn phương trình Poisson:

2φ = 4π2Gρ, 

trong đó ρ(r) là mật độ khối lượng, còn G là hằng số trọng trường.

Bài toán 1 chiều với điều kiện biên Dirichlet

Ta xét nghiệm của phương trình Poisson trong không gian một chiều làm ví dụ đơn giản. Giả sử rằng

d2u(x) / dx2 = v(x), 

với xl ≤ x ≤ xh, tuân theo điều kiện biên Dirichlet u(xl) = ulu(xh) = uh.

Bước đầu tiên là chia miền xl ≤ x ≤ xh thành những đoạn dài bằng nhau sao cho các đỉnh nằm ở những điểm nút lưới

xi = xl + i (xh − xl) / (N + 1),    [e57]

với i = 1, N. Các biên, xlxh, lần lượt tương ứng với i = 0i = N + 1.

Tiếp theo, ta làm rời rạc biểu thức vi phân d2u / dx2 trên lưới. Cách làm rời rạc dễ thấy nhất là

d2u(xi) / dx2 = (ui − 1 − 2 ui + ui + 1)/(δx)2  +  O(δx)2.    [e58]

Ở đây, δx = (xh − xl) / (N + 1), và ui ≡ u(xi). Kiểu làm rời rạc này được gọi là sơ đồ sai phân trung tâm bậc hai. Nó được gọi là “bậc hai” vì sai số chặt cụt là O(δx)2, như có thể dễ dàng cho thấy từ khai triển chuỗi Taylor. Dĩ nhiên, một sơ đồ bậc n sẽ có sai số chặt cụt là O(δx)n. Đó là một sơ đồ “sai phân trung tâm” vì nó đối xứng quanh nút trung tâm, xi. Dạng phương trình Poisson sau khi làm rời rạc là

ui − 1 − 2 ui + ui + 1 = vi (δx)2, 

với i = 1, N, trong đó vi ≡ v(xi). Hơn nữa, u0 = uluN + 1 = uh.

Sẽ rất có ích nếu ta coi hệ phương trình sau khi rời rạc trên dưới dạng phương trình ma trận. Đặt u = (u1,  u2,  ⋯,  uN) là véc-tơ chứa các giá trị của u, và đặt

w = [v1 (δx)2 − ul,  v2 (δx)2,  v3 (δx)2,  ⋯,  vN − 1 (δx)2,  vN (δx)2 − uh]   [e510]

là véc-tơ chứa các số hạng nguồn. Các phương trình rời rạc có thể được viết thành

Mu = w.     [e511]

Ma trận M có dạng

\displaystyle{ \textbf{M} = \left(\begin{array}{rrrrrr}  -2 & 1 & 0 & 0 & 0 & 0 \\  1 & -2 & 1 & 0 & 0 & 0 \\  0 & 1 & -2 & 1 & 0 & 0 \\  0 & 0 & 1 & -2 & 1 & 0 \\  0 & 0 & 0 & 1 & -2 & 1 \\  0 & 0 & 0 & 0 & 1 & -2 \\  \end{array}  \right) }

với N = 6. Dễ thấy ta có thể khái quát đối với những giá trị khác của N. Ma trận M được gọi là ma trận ba đường chéo, vì ở đó chỉ có những phần tử nằm trên ba đường chéo chính mới khác không.

Nghiệm hình thức của PT ([e511]) là

u = M − 1w, 

trong đó M − 1 là ma trận nghịch đảo của M. Không may là thuật toán tổng quát có hiệu quả nhất để tính nghịch đảo của một ma trận N × N—phương pháp khử Gauss-Jordan với chốt từng phần—vẫn cần đến O(N3) phép tính số học. Dễ thấy rằng thang bậc tỉ lệ này thật tai hại khi phải tìm nghiệm sai phân hữu hạn của phương trình Poisson. Mỗi khi ta tăng gấp đôi độ phân giải (nghĩa là gấp đôi số điểm nút) thì thời gian chạy CPU tăng gấp tám lần. Hệ quả là, việc thêm vào một chiều thứ hai (vốn yêu cầu số điểm nút phải được bình phương lên) thí sẽ cực kì hao tốn thời gian chạy CPU đến mức không chấp nhận được. Tuy nhiên, có một cách làm quen thuộc để lấy nghịch đảo một ma trận ba đường chéo N × N mà chỉ cần O(N) phép tính số học.

Xét phương trình cho một ma trận ba đường chéo N × N tổng quát: Mu = w. Đặt a, b, và c lần lượt là các véc-tơ phía trái, giữa, và bên phải của đường chéo chính trong ma trận. Lưu ý rằng a1cN đều không xác định, và để cho tiện ta có thể cho chúng bằng không. Bây giờ phương trình ma trận có thể viết dưới dạng

aiui − 1 + biui + ciui + 1 = wi,    [e515]

với i = 1, N. Hãy tìm một nghiệm có dạng

ui + 1 = xiui + yi.    [e517]

Thay thế vào PT ([e515]) cho ta

aiui − 1 + biui + ci (xiui + yi) = wi, 

vốn có thể được sắp xếp để cho

\displaystyle{u_i = - \frac{a_i\,u_{i-1}}{b_i+c_i\,x_i} + \frac{w_i-c_i\,y_i}{b_i+c_i\,x_i}.}

Tuy nhiên, nếu PT ([e517]) là dạng chung thì ta có thể viết ui = xi − 1ui − 1 + yi − 1. So sánh với phương trình trước ta được

\displaystyle{ x_{i-1} = - \frac{a_i}{b_i+c_i\,x_i},}     [e520]

\displaystyle{y_{i-1} = \frac{w_i-c_i\,y_i}{b_i+c_i\,x_i}.}     [e521]

Bây giờ ta có thể giải phương trình ma trận ba đường chéo theo hai bước. Ở bước thứ nhất, ta dò ngược theo đường chéo trước từ i = N đến 1 bằng các PT ([e520]) và ([e521]). Như vậy,

\displaystyle{x_{N-1} = -\frac{a_N}{b_N},\mbox{\hspace{1cm}}y_{N-1} = \frac{w_N}{b_N},}

cN = 0. Hơn nữa,

\displaystyle{ x_i = - \frac{a_{i+1}}{b_{i+1}+c_{i+1}\,x_{i+1}},\mbox{\hspace{1cm}}y_i =  \frac{w_{i+1}-c_{i+1}\,y_{i+1}}{b_{i+1}+c_{i+1}\,x_{i+1}} }

với i = N − 2, 1. Cuối cùng,

\displaystyle{ x_0 = 0,\mbox{\hspace{1cm}}y_0 = \frac{w_1-c_1\,y_1}{b_1+c_1\,x_1}, }

a1 = 0. Bây giờ ta đã định nghĩa được tất cả xiyi. Ở bước thứ hai, ta dò xuôi theo đường chéo trước từ i = 0 đến N − 1 bằng PT ([e517]). Như vậy,

u1 = y0, 

x0 = 0, và

ui + 1 = xiui + yi

với i = 1, N − 1. Ta đã nghịch đảo được ma trận ba đường chéo để giải phương trình chỉ bằng O(N) phép tính số học.

Rõ ràng là ta có thể dùng thuật toán trên để nghịch đảo PT ([e511]), với các số hạng nguồn cho trong PT ([e510]), và các đường chéo của ma trận M cho bởi ai = 1 với i = 2, N, cùng bi =  − 2 với i = 1, N, và ci = 1 với i = 1, N − 1.

Ví dụ về đoạn chương trình giải ma trận ba đường chéo

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

// Tridiagonal.cpp 

// Function to invert tridiagonal matrix equation. 
// Left, centre, and right diagonal elements of matrix 
// stored in arrays a, b, c, respectively. 
// Right-hand side stored in array w. 
// Solution written to array u. 

// Matrix is NxN. Arrays a, b, c, w, u assumed to be of extent N+2, 
// with redundant 0 and N+1 elements. 

#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) 
{ 
  // Find N. Declare local arrays.     
  int N = a.extent(0) - 2; 
  Array<double,1> x(N), y(N); 

  // Scan up diagonal from i = N to 1 
  x(N-1) = - a(N) / b(N); 
  y(N-1) = w(N) / b(N); 
  for (int i = N-2; i > 0; i--) 
    { 
      x(i) = - a(i+1) / (b(i+1) + c(i+1) * x(i+1)); 
      y(i) = (w(i+1) - c(i+1) * y(i+1)) / (b(i+1) + c(i+1) * x(i+1)); 
    } 
  x(0) = 0.; 
  y(0) = (w(1) - c(1) * y(1)) / (b(1) + c(1) * x(1)); 

  // Scan down diagonal from i = 0 to N-1 
  u(1) = y(0); 
  for (int i = 1; i < N; i++) 
    u(i+1) = x(i) * u(i) + y(i); 
}

Bài toán 1 chiều với điều kiện biên hỗn hợp

Trước đây, ta đã giải phương trình Poisson trong không gian một chiều với điều kiện biên Dirichlet, vốn là điều kiện biên đơn giản nhất có thể. Bây giờ, hãy xét điều kiện biên tổng quát hơn nhiều như sau:

αlu(x) + βldu(x) / dx = γl,     [e527]

tại x = xl, và

αhu(x) + βhdu(x) / dx = γh,     [e528]

tại x = xh. Ở đây, αl, βl, v.v., là các hằng số đã biết. Điều kiện biên trên được gọi là hỗn hợp, vì chúng là sự kết hợp giữa các điều kiện biên Dirichlet và Neumann.

Dùng kí hiệu như trên, có thể viết dạng rời rạc của các PT ([e527]) và ([e528]) lần lượt như sau:

\displaystyle{ \begin{array} {ccc}  \alpha_l\,u_0 + \beta_l\,\frac{u_1-u_0}{\delta x} &=& \gamma_l,\\[0.5ex]  \alpha_h\,u_{N+1} + \beta_h\,\frac{u_{N+1}-u_N}{\delta x} &=& \gamma_h.  \end{array} }      [e529] và [e530]

Có thể sắp xếp lại các biểu thức trên để thu được

\displaystyle{ u_0 = \frac{\gamma_l\,\delta x-\beta_l\,u_1}{\alpha_l\,\delta x -\beta_l},}    [e531]

\displaystyle{ u_{N+1} = \frac{\gamma_h\,\delta x + \beta_h\,u_N}{\alpha_h\,\delta x +\beta_h}.}    [e532]

Dùng các PT ([e58]), ([e531]), và ([e532]), có thể giản hóa bài toán về một phương trình ma trận ba đường chéo Mu = w, trong đó các đường chéo trái, giữa và phải của M chứa những phần tử lần lượt là ai = 1 với i = 2, N,

\displaystyle{ b_1 = -2 - \frac{\beta_l}{\alpha_l\,\delta x-\beta_l},}

bi =  − 2 với i = 2, N − 1,

\displaystyle{ b_N = -2 + \frac{\beta_h}{\alpha_h\,\delta x+\beta_h},}

ci = 1 với i = 1, N − 1. Các phần tử của vế phải là

\displaystyle{ w_1 = v_1\,(\delta x)^2 - \frac{\gamma_l\,\delta x}{\alpha_l\,\delta x-\beta_l}, }

wi = vi (δx)2 với i = 2, N − 1, và

\displaystyle{ w_N = v_N\,(\delta x)^2 - \frac{\gamma_h\,\delta x}{\alpha_h\,\delta x+\beta_h}. }

Có thể nghịch đảo phương trình ma trận ba đường chéo này bằng thuật toán đã trình bày ở trên.

Ví dụ về đoạn chương trình giải bài toán Poisson một chiều

Dưới đây là một đoạn chương trình giải bài toán Poisson một chiều trong đó có dùng cách giải ma trận ba đường chéo từ trước và thư viện Blitz++ (xem Mục [mảng nhiều chiều, ch2]).

// Poisson1D.cpp 

// Function to solve Poisson's equation in 1-d: 

//  d^2 u / dx^2 = v  for  xl <= x <= xh 

//  alpha_l u + beta_l du/dx = gamma_l  at x=xl 

//  alpha_h u + beta_h du/dx = gamma_h  at x=xh 

// Arrays u and v assumed to be of extent N+2. 

// Now, ith elements of arrays correspond to 

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

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

#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 Poisson1D (Array<double,1>& u, Array<double,1> v, 
                double alpha_l, double beta_l, double gamma_l, 
                double alpha_h, double beta_h, double gamma_h, 
                double dx) 
{ 
  // 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) = 1.; 
  for (int i = 1; i <= N; i++) b(i) = -2.; 
  b(1) -= beta_l / (alpha_l * dx - beta_l); 
  b(N) += beta_h / (alpha_h * dx + beta_h); 
  for (int i = 1; i <= N-1; i++) c(i) = 1.; 

  // Initialize right-hand side vector 
  for (int i = 1; i <= N; i++) 
      w(i) = v(i) * dx * dx; 
  w(1) -= gamma_l * dx / (alpha_l * dx - beta_l); 
  w(N) -= gamma_h * dx / (alpha_h * dx + beta_h); 

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

  // Calculate i=0 and i=N+1 values 
  u(0) = (gamma_l * dx - beta_l * u(1)) / 
    (alpha_l * dx - beta_l); 
  u(N+1) = (gamma_h * dx + beta_h * u(N)) / 
    (alpha_h * dx + beta_h); 
}

Ví dụ về kết quả của phương trình Poisson một chiều

Bây giờ, ta hãy giải phương trình Poisson trong một chiều, với các điều kiện biên hỗn hợp, bằng kĩ thuật sai phân hữu hạn như đã trình bày. Ta cần tìm nghiệm dưới dạng

d2u(x) / dx2 = v(x), 

trong miền 0 ≤ x ≤ 1, với v(x) = 1 − 2 x2. Các điều kiện biên tại xl = 0xh = 1 có dạng hỗn hợp như ở các PT ([e527]) và ([e528]). Dĩ nhiên là ta có thể giải bài toán này theo cách giải tích. Thật ra,

u(x) = g + hx + x2 / 2 − x4 / 6, 

trong đó

\displaystyle{ g = \frac{\gamma_l\,(\alpha_h+\beta_h) - \beta_l\,[\gamma_h-(\alpha_h+\beta_h)/3]}  {\alpha_l\,\alpha_h + \alpha_l\,\beta_h-\beta_l\,\alpha_h},}

\displaystyle{ h = \frac{ \alpha_l\,[\gamma_h-(\alpha_h+\beta_h)/3]-\gamma_l\,\alpha_h}  {\alpha_l\,\alpha_h + \alpha_l\,\beta_h-\beta_l\,\alpha_h}.}

Hình [f5x1] cho thấy sự so sánh giữa các nghiệm giải tích và sai phân hữu hạn cho trường hợp N = 100. Có thể thấy rằng nghiệm sai phân hữu hạn gần như phản ánh đúng nghiệm giải tích.

Hình [f5x1]. Nghiệm của phương trình Poisson trong không gian một chiều với v = 1 − 2 x2, αl = 1, βl =  − 1, γl = 1, αh = 1, βh = 1, và γh = 1. Đường nét đứt (bị che khuất) biểu thị nghiệm giải tích, còn những hình tam giác biểu thị nghiệm sai phân hữu hạn cho trường hợp N = 100

Bài toán 2 chiều với điều kiện biên Dirichlet

Ta hãy xét lời giải bài toán Poisson trong không gian hai chiều. Giả sử rằng

2u(x, y)/∂x2 + ∂2u(x, y)/∂y2 = v(x, y),  [e541]

với xl ≤ x ≤ xh, và 0 ≤ y ≤ L. Cũng tương tự như cách dùng cho bài toán một chiều trước đây, ta có thể làm rời rạc phương trình hai chiều trên bằng một lược đồ sai phân trung tâm bậc hai, theo cả hai chiều xy. Thật không may là lược đồ này sẽ cho ra một hệ phương trình mà ta không thể giản hóa về một phương trình ma trận ba đường chéo. Thực ra, tất cả các thuật toán số trị hiệu quả để giải dạng bài toán này đều có tính chất lặp. Chẳng hạn, phương pháp Jacobi, phương pháp Gauss-Seidel, phương pháp giãn quá liên tiếp, và phương pháp nhiều lưới.2 Thật tiếc là những phương pháp đó, nếu không phải cực kì phức tạp như phương pháp nhiều lưới (và do đó vượt ngoài phạm vi khóa học này), thì lại hội tụ rất kém. Trong phần tiếp theo, thay vì tìm hiểu những phương pháp lặp không được hiệu quả, ta sẽ tìm hiểu về một phương pháp không lặp nhưng hoạt động hiệu quả với một nhóm các bài toán cụ thể. Phương pháp này được gọi là phương pháp phổ, vì nó liên quan đến việc khai triển uv dưới dạng chuỗi Fourier bị chặt cụt theo phương y.

Giả sử rằng u(x, y) thỏa mãn điều kiện biên hỗn hợp theo phương x, nghĩa là

αlu(x, y) + βl ∂u(x, y) / ∂x = γl(y),   [e528a]

tại x = xl, và

αhu(x, y) + βh ∂u(x, y) / ∂x = γh(y),   [e528a]

tại x = xh. Ở đây, αl, βl, v.v., là các hằng số đã biết, còn γl, γh là các hàm theo y đã biết. Hơn nữa, giả sử rằng u(x, y) thỏa mãn điều kiện biên Dirichlet đơn giản sau theo phương y:

u(x, 0) = u(x, L) = 0.      [e5.44]

Lưu ý rằng, vì u(x, y) là một thế, và do đó có tính bất định [ta có thể tùy ý cộng hoặc trừ một hằng số], điều kiện biên trên tương đương với việc yêu cầu u phải nhận cùng giá trị hằng số ở cả biên trên và biên dưới theo phương y.

Ta hãy viết u(x, y) dưới dạng chuỗi Fourier theo phương y:

u(x, y) = ∑ j = 0Uj(x) sin(jπy / L).        [e545]

Lưu ý rằng biểu thức trên của u tự thỏa mãn các điều kiện biên theo phương y. Các hàm sin(jπy / L) có tính trực giao, và hình thành một bộ đầy đủ trong khoảng y = 0, L. Thật ra,

(2 / L)∫ 0Lsin(jπy / L) sin(kπy / L) dy = δjk. 

Vì vậy, ta có thể viết số hạng nguồn dưới dạng

v(x, y) = ∑ j = 0Vj(x) sin(jπy / L),       [e547]

trong đó

Vj(x) = (2 / L)∫ 0Lv(x, y) sin(jπy / L) dy.       [ffta1]

Hơn nữa, các điều kiện biên theo phương x trở thành

αlUj(x) + βldUj(x) / dx = Γ lj,      [e549]

tại x = xl, và

αhUj(x) + βhdUj(x) / dx = Γ hj,       [e550]

tại x = xh, trong đó

Γ lj = (2 / L)∫ 0Lγl(y) sin(jπy / L) dy,      [ffta2]

v.v.  Thay các PT ([e545]) và ([e547]) vào PT ([e541]), rồi cân bằng những hệ số của sin(jπy / L) (vì các hàm này trực giao nhau), ta có

\displaystyle{ \frac{d^2 U_j(x)}{dx^2} - \frac{j^2\,\pi^2}{L^2}\,U_j(x) = V_j(x), }

với j = 0, ∞. Bây giờ ta có thể làm rời rạc phương trình theo phương y bằng cách chặt cụt biểu thức khai triển Fourier này, nghĩa là bằng cách chỉ giải các phương trình trên cho j = 0, J, thay vì j = 0, ∞. Điều này hoàn toàn tương đương với việc làm rời rạc theo phương y trên lưới có các nút cách đều nhau yj = jL / J. Phương trình được làm rời rạc theo phương x bằng cách chia miền thành những đoạn bằng nhau, theo PT ([e57]), rồi xấp xỉ d2 / dx2 theo lược đồ sai phân trung tâm bậc hai như ở PT ([e58]). Theo đó, ta có

Ui − 1, j − (2 + j2κ2) Ui, j + Ui + 1, j = Vi, j (δx)2,       [e5a]

với i = 1, N and j = 0, J. Ở đây, Ui, j ≡ Uj(xi), Vi, j ≡ Vj(xi), và κ = πδx / L. Các điều kiện biên ([e549]) và ([e550]) được rời rạc để cho

\displaystyle{ U_{0,j} = \frac{\Gamma_{l\,j}\,\delta x-\beta_l\,U_{1,j}}{\alpha_l\,\delta x -\beta_l},}      [e5b]

\displaystyle{ U_{N+1,j} = \frac{\Gamma_{h\,j}\,\delta x + \beta_h\,U_{N,j}}{\alpha_h\,\delta x +\beta_h},}    [e5c]

với j = 0, J. Các PT ([e5a]), ([e5b]), và ([e5c]) lập thành J + 1 phương trình ma trận tam giác độc lập nhau (mỗi phương trình cho một giá trị j riêng). Các phương trình này có thể được nghịch đảo bằng thuật toán đã trình bày ở Mục [điều kiện biên phức hợp], để cho ra các Ui, j. Cuối cùng, những giá trị u(xi, yj) có thể được tái lập từ PT ([e545]). Như vậy là ta đã giải được bài toán.

Bài toán 2 chiều với điều kiện biên Neumann

Ta hãy thực hiện lại những tính toán trên nhưng thay điều kiện biên Dirichlet ([e5.44]) với điều kiện biên Neumann đơn giản sau:

\displaystyle{ \frac{\partial u(x,y=0)}{\partial y} = \frac{\partial u(x,y=L)}{\partial y} = 0. }

Trong trường hợp này, có thể biểu diễn u(x, y) dưới dạng

u(x, y) = ∑ j = 0Uj(x) cos(jπy / L),     [e5x]

vốn tự thỏa mãn các điều kiện biên theo phương y. Tương tự, ta có thể viết số hạng nguồn v(x, y) theo dạng

v(x, y) = ∑ j = 0Vj(x) cos(jπy / L), 

trong đó

Vj(x) = (2 / L)∫ 0Lv(x, y) cos(jπy / L) dy,      [nb1]

(2 / L)∫ 0Lcos(jπy / L) cos(kπy / L) dy = δjk. 

Cuối cùng, các điều kiện biên theo phương x trở thành

αlUj(x) + βldUj(x) / dx = Γ lj, 

tại x = xl, và

αhUj(x) + βhdUj(x) / dx = Γ hj, 

tại x = xh, trong đó

Γ lj = (2 / L)∫ 0Lγl(y) cos(jπy / L) dy,      [nb2]

v.v. Tuy nhiên, cần lưu ý rằng nhân tử đứng trước các tích phân trong PT ([nb1]) và ([nb2]) nhận giá trị đặc biệt 1 / L với thành phần điều hòa j = 0. Cũng như trước, ta chặt cụt chuỗi khai triển Fourier theo phương y, rồi làm rời rạc theo phương x, để thu được hệ phương trình ma trận ba đường chéo như trong các PT ([e5a]), ([e5b]), và ([e5c]). Có thể giải được các phương trình này để thu được Ui, j, rồi sau đó tái lập các u(xi, yj) từ PT ([e5x]). Như vậy, ta đã giải được bài toán.


  1. Ở tài liệu gốc, simply-connected là không gian tô-pô trong đó hai điểm bất kì trong không gian được nối bằng đường nằm trọn trong chính không gian đó.
  2. Xem 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), Sect. 19.5.

18 phản hồi

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

18 responses to “Chương 5: Phương trình Poisson (Phần 1)

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

  2. Pingback: Chương 5: Phương trình Poisson (Phần 2) | Blog của Chiến

  3. Joke

    Cảm ơn bạn về bài viết. Mình cũng đang giải phương trình dạng này nhưng trong môi trường nước.

    • Chúc mừng bạn! Hi vọng rằng từ bài viết bạn sẽ tìm được những thông tin có ích để lập trình giải được bài toán thuỷ động.

  4. Joke

    Bạn Quang Chiến vui lòng cho mình hỏi chút. Trong một bài toán giải phương trình khuếch tán, một biên nào đó (giả sử biên số 1 chẳng hạn) nằm trên mặt bao S của phần thể tích V mình đang xét, liệu có thể đồng thời áp dụng cả điều kiện biên Neumann và điều kiện biên Dirichlet cho biên đó được không vậy bạn (xét theo mặt toán – lý, chứ không nói về việc lập trình)?

  5. Về nguyên tắc là không được vì như thế sẽ khiến cho số phương trình dư ra so với số ẩn (thừa 1 phương trình) khiến cho hệ không xác định. Về vật lý, chảng hạn xét phương trình truyền nhiệt, thì nếu ta ấn định nhiệt độ ở mép một tấm kim loại phẳng thì tự nó sẽ quy định thông lượng nhiệt truyền vào vùng trong, tức là quy định luôn giá trị gra-đi-en nhiệt rồi.

    • Không nên nhầm với việc dùng điều kiện biên Robin có cả giá trị và đạo hàm, vì thực ra điều kiện biên này cũng chỉ có 1 phương trình (chứ không phải là hai phương trình Dirichlet và Neumann).

      • Joke

        Theo mình hiểu, thì điều kiện biên Robin là sự kết hợp (theo toán tử cộng, contribute) cả điều kiện biên Neumann và Dirichlet vào một biên nào đó.

        Tuy nhiên nếu xét theo ví dụ của bạn Chiến, tại một mép nào đó của tấm kim loại (giả sử gọi đó là biên 1) ta quy định thông lượng của nó bằng 0 (Neumann) và sau đó cũng tại biên 1 đấy, ta quy định nhiệt độ cố định bằng 50oC chẳng hạn, thì bài toán hoàn toàn có thể hợp lý (vì nhiệt sẽ truyền tuân theo định luật khuếch tán Fick II, vì mình áp dụng định định luật này trong toàn bộ thể tích V đang xét, thì thể tích V sẽ phải bao gồm biên 1 chứ ạ)?
        Nếu mình có điều gì sai, mong được chỉ bảo ạ.
        Cảm ơn anh Chiến.

  6. Joke

    Mình đang giải bài toán này bằng dùng phần mềm mô phỏng Comsol, và bài thắc mắc của mình, cũng được hỏi trên diễn đàn của phần mềm này cách đây nửa tiếng, bạn Chiến vui lòng xem qua cho hiểu thêm ý câu hỏi của mình. Cảm ơn bạn nhiều.
    http://www.comsol.com/community/forums/general/thread/31800/#p86840

    • Cám ơn bạn đã đưa thông tin cụ thể. Dù vẫn chưa có được câu trả lời thỏa đáng nhưng hiện giờ mình thấy thế này:
      * Cái thể tích V chỉ chứa các điểm bên trong miền tính toán thôi, những điểm giới hạn bởi (1,1), (1,N), (N,1), (N,N). Còn biên DIrichlet thì áp dụng cho những điểm nút ảo (“ghost”) chẳng hạn từ (0,0)–(0,N+1)–(N+1,N+1)–(N+1,0). Biên Neumann ứng với luồng nhiệt truyền từ điểm ghost vào trong, hoặc ngược lại; chẳng hạn như giữa hai điểm (3,N) và (3,N+1).
      * Thực ra thì mình vừa xem thấy biên Cauchy gồm cả hai phương trình Dirichlet và Neumann, áp dụng được vào bài toán truyền nhiệt. Nhưng lại để ý thấy câu hỏi ban đầu của bạn là phương trình khuếch tán (PT dạng hyperbol); nó khác với phương trình truyền nhiệt (parabol). Liệu đây có phải là vấn đề không? Trên diễn đàn Comsol bạn cần nói cụ thể hơn là đang giải PT loại nào.
      * Riêng về bài toán truyền nhiệt có vẻ giả thiết bạn đưa ra là hợp lý. Mình xét trường hợp đơn giản là thanh có 2 đầu gắn vào vách ở các nhiệt độ khác nhau là TA và TB (Dirichlet). Vậy hiển nhiên sau một thời gian phân bố nhiệt trên thanh là dạng đường thẳng nối các giá trị TA và TB rồi. Luồng nhiệt chạy đều trên toàn thanh, với giá trị bằng độ dốc của đường phân bố này. Nếu bạn vẫn muốn đưa vào gra-đi-en nhiệt tại B thì nó sẽ khiến cho phân bố nhiệt trên thanh thay đổi; nói cách khác cả giá trị TB và (dT/dx)B đều có ảnh hưởng đến kết quả. Bạn có thể xem: http://50.57.233.96/comput_phys/heat1.png

  7. phamquyen

    ban ah.ban co phan nhung dieu kien vat ly dan den phuong trinh poisson ko

    • Mình chưa hiểu rõ ý bạn lắm. Còn nếu bạn hỏi về ý nghĩa vật lý của phương trình Poisson thì, như đoạn đầu bài viết này, điện trường và trọng trường có thể được biểu diễn bởi phương trình Poisson. Tổng quát hơn, một trường thế có thể được biểu diễn bởi phương trình Poisson. Trong cơ học thủy khí, động thái của một chất lưu (fluid) lý tưởng cũng có thể biểu diễn bởi phương trình Poisson; ở đây \varphi là áp suất p.

      • phamquyen

        minh dang lam de tai lien quan den phuong trinh poisson va trong phan co ban cua minh can phai lam “Những điều kiện vật lý dẫn đến phương trình poisson” .Neu ban co tai lieu ve phan nay cho minh xin nhe.Cam on ban

      • Trong số những quyển sách tựa đề kiểu như (mathematical, method, physics) thì không đâu nói kĩ về điều kiện vật lý của PT Poisson ? Ví dụ cốt yếu của họ chỉ là về điện trường tĩnh như vậy thôi. Các ứng dụng khác về nhiệt, cơ lỏng, v.v. thì họ cứ nghiễm nhiên áp dụng luôn chẳng đề cập gì đến “điều kiện” cả. Đây là kiến thức cơ bản về lý thuyết trường rồi.

  8. bạn phạm quyền ơi, bạn đã tìm được tài liệu về những điều kiện vật lý dẫn đến phương trình poisson chưa bạn? bạn có thể chia sẻ cho mình được ko? mình cũng đang làm đề tài tốt nghiệp thì có 1 phần liên quan tới điều kiện đó, bạn có thì cho mình xin với nhé. tks bạn

  9. Bùi Xuân Trường

    Bạn Quangchien cho mình hỏi chút là phương trình Poisson cho bài toán phân bố thế của hệ thống nối đất HVDC có dạng như nào vậy bạn?

  10. Bùi Xuân Trường

    Cảm ơn bạn nhiều nhé, mình cũng tìm hiểu rồi mà chưa thấy, có thông tin gì mới thì cho mình biết nhé.

Gửi phản hồ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