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

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

Phương trình sóng 1 chiều

Xét một sóng điện từ phân cực phẳng truyền trong chân không dọc theo trục z. Coi rằng các trường điện và từ có dạng E = [Ex(z, t), 0, 0], và B = [0, By(z, t), 0]. Bây giờ, theo định luật Maxwell,

\begin{array}{ccc}  \partial E_x / \partial t + c\, \partial H_y / \partial z &=& 0,\\  \partial H_y / \partial t + c\, \partial E_x / \partial z &=& 0,\\  \end{array}             [wave1da] & [wave1db]

trong đó Hy = cBy, và c là vận tốc ánh sáng. Lưu ý rằng các phương trình trên có dạng hệ hai phương trình chuyển tải. Ta hãy đi tìm nghiệm số trị của các phương trình này trong một vùng 0 ≤ z ≤ L nào đó được bao kín bởi những vách cách điện tại z = 0z = L. Giờ đây, cả điện trường tiếp tuyến và từ trường pháp tuyến phải đều bằng không tại một vật dẫn tuyệt đối. Từ đó ta có

Ex(0, t) = Ex(L, t) = 0.              [wavebca]

Bên cạnh đó, PT ([wave1da]) cho ta

Hy(0, t) / ∂z = ∂Hy(L, t) / ∂z = 0.                  [wavebcb]

Điều ta trông đợi ở các PT ([wave1da]) và ([wave1db]) là khả năng tính được các nghiệm dạng sóng mà nảy qua lại giữa hai vách cách điện.

Cũng như trước, 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, ⋯. Mặt khác, theo phương z, ta làm rời rạc trên luới đều zi = iδz, với i = 0, I, trong đó δz = L / I. Dùng một lược đồ sai phân Crank-Nicolson theo thời gian giống với lược đồ đề cập ở Mục [cnwave], các PT ([wave1da])–([wave1db]) cho ta

\begin{array}{ccc}  \displaystyle{ \frac{(E_x)_i^{n+1}-(E_x)_i^n}{\delta t} +\frac{c}{2}\left(\frac{\partial H_y}{\partial z}\right)_i^{n+1}  + \frac{c}{2}\left(\frac{\partial H_y}{\partial z}\right)_i^n } &=& 0, \\  \displaystyle{ \frac{(H_y)_i^{n+1}-(H_y)_i^n}{\delta t} +\frac{c}{2}\left(\frac{\partial E_x}{\partial z}\right)_i^{n+1}  + \frac{c}{2}\left(\frac{\partial E_x}{\partial z}\right)_i^n } &=& 0,\\  \end{array}           [cnwave1da] & [cnwave1db]

trong đó (Ex)in ≡ Ex(zi, tn), v.v.

Dùng phương pháp Fourier, ta viết

\begin{array}{ccc}  (E_x)_i^n &=& \sum_{j=0,I} \hat{E}_j^n\,\sin(i\,j\,\pi/I),\\  (H_y)_i^n &=& \sum_{j=0,I} \hat{H}_j^n\,\cos(i\,j\,\pi/I),\\  \end{array}

vốn tự thỏa mãn các điều kiện biên ([wavebca]) và ([wavebcb]). Các PT ([cnwave1da]) và ([cnwave1db]) cho ra

\begin{array}{ccc}  \hat{E}_i^{n+1}-\hat{E}_i^n - i\,D\,(\hat{H}_i^{n+1}+\hat{H}_i^n) &=& 0,\\  \hat{H}_i^{n+1}-\hat{H}_i^n + i\,D\,(\hat{E}_i^{n+1}+\hat{E}_i^n) &=& 0,\\  \end{array}

với D = πcδt / (2 L). Từ đó dẫn đến

\displaystyle{  \hat{E}_i^{n+1} = +\frac{2\,i\,D}{1+i^2\,D^2}\,\hat{H}_i^n + \frac{1-i^2\,D^2}{1+i^2\,D^2}\,\hat{E}_i^n,}

\displaystyle{  \hat{H}_i^{n+1} = -\frac{2\,i\,D}{1+i^2\,D^2}\,\hat{E}_i^n + \frac{1-i^2\,D^2}{1+i^2\,D^2}\,\hat{H}_i^n,}                    [stepa] & [stepb]

với i = 0, I.

Ta hãy tìm giá trị riêng λ của hệ tuyến tính trên, với giả thiết rằng (\hat{E}_{i}^{n+1},\hat{H}_{i}^{n+1}) = \lambda(\hat{E}_{i}^{n},\hat{H}_{i}^{n}). Sau khi tính toán, ta được

\displaystyle{  \lambda^2 - 2\,\frac{1-i^2\,D^2}{1+i^2\,D^2}\,\lambda + 1 = 0}.

Với mọi i, hai nghiệm của phương trình bậc hai trên đều là số phức, nhưng có độ lớn đều bằng đơn vị, hay λ∣ = 1. Điều này nghĩa là lược đồ sai phân hiện giờ có cả tính ổn định số trị (nếu λ∣ > 1 thì lược đồ sẽ bất ổn định, vì hàm riêng tương ứng sẽ bị khuếch đại sau mỗi bước thời gian) và tính không phân tán (nếu λ∣ < 1 thì tất cả hàm điều hòa Fourier cuối cùng sẽ tiêu biến, và do đó nghiệm cũng vậy). Lưu ý rằng phép tính trên tương đương với phân tích ổn định von Neumann.

Đoạn chương trình dưới đây giải phương trình sóng 1 chiều bằng lược đồ Crank-Nicolson đã đề cập ở trên. Chương trình trước hết biến đổi Fourier các đại lượng ExHy, chọn bước thời gian dựa vào các PT ([stepa]) và ([stepb]), rồi tái lập lại ExHy qua phép biến đổi Fourier ngược.

// Wave1D.cpp 

// Function to evolve 1-d wave equation: 

//  d E_x / dt + c d H_y / dz = 0 

//  d H_y / dt + c d E_x / dz = 0 

// in region 0 &lt; z &lt; L.  

// Boundary conditions:  

//  E_x(0,t) = E_x(L,t) = 0 

//  d H_y(0,t) / dz =  d H_y(L,t) / dz = 0 

// Arrays Ex, Hy assumed to be of extent I+1. 

// Now, ith elements of arrays correspond to 

//  z_i = i * L / J      i=0,I 

// Also, D = pi c dt / (2 L), where dt is time-step. 

// Uses Crank-Nicolson scheme. 

#include &lt;blitz/array.h&gt; 

using namespace blitz; 

void fft_forward_cos (Array&lt;double,1&gt; f, Array&lt;double,1&gt;&amp; F); 
void fft_backward_cos (Array&lt;double,1&gt; F, Array&lt;double,1&gt;&amp; f); 
void fft_forward_sin (Array&lt;double,1&gt; f, Array&lt;double,1&gt;&amp; F); 
void fft_backward_sin (Array&lt;double,1&gt; F, Array&lt;double,1&gt;&amp; f); 

void Wave1D (Array&lt;double,1&gt;&amp; Ex, Array&lt;double,1&gt;&amp; Hy, double D) 
{ 
  // Find I. Declare local arrays 
  int I = Ex.extent(0) - 1; 
  Array&lt;double,1&gt; EE(I+1), HH(I+1), EE0(I+1), HH0(I+1); 

  // Fourier transform Ex and Hy 
  fft_forward_sin (Ex, EE0); 
  fft_forward_cos (Hy, HH0); 

  // Evolve EE and HH 
  for (int i = 0; i &lt;= I; i++) 
    { 
      double x = double (i) * D; 
      double fp = 1. + x*x; 
      double fm = 1. - x*x; 

      EE(i) = 2. * x * HH0(i) + fm * EE0(i); 
      EE(i) /= fp; 
      HH(i) = - 2. * x * EE0(i) + fm * HH0(i); 
      HH(i) /= fp; 
    } 

  // Reconstruct Ex and Hy via inverse Fourier transform 
  fft_backward_sin (EE, Ex); 
  fft_backward_cos (HH, Hy); 
}

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

Ex(z, 0) = Hy(z, 0) = exp[ − 100 (z − 0. 5)2], 

và tính toán được thực hiện với L = 1, c = 1, δt = 5 × 10−3, và N = 100. Có thể thấy rằng xung chuyển động sang bên phải, phản xạ ở vách cách điện đặt ở z = 1, rồi chuyển động về phía trái. Lưu ý rằng Ex =  +Hy khi sóng còn xa vách cách điện và đang chuyển động qua phải, còn Ex =  −Hy khi sóng còn xa vách cách điện và đang chuyển động qua trái.

wave1d

Hình. Sự lan truyền của xung sóng Gauss 1 chiều. Tính toán số trị được thực hiện với c = 1, δt = 5 × 10−3, và N = 100. Đường liền nét biểu thị Ex, còn đường gạch đứt biểu thị Hy. Các biểu đồ con ở trên-trái, trên-phải, ngang-trái, ngang-phải, dưới-trái, và dưới-phải lần lượt biểu thị nghiệm lúc t = 0. 0, 0. 2, 0. 4, 0. 6, 0. 8, và 1. 0. [fwave1d]

Hình [wave1da] cho thấy một ví dụ tính toán thứ hai được thực hiện với các tham số giống như ví dụ đầu. Trong tính toán lần này, xung được phép phản xạ ở vách cách điện mười lần trước khi trở về vị trí ban đầu của nó. Lưu ý rằng biên độ và hình dáng của xung gần như vẫn giữ nguyên suốt quá trình này. Việc xung gần như trở về đúng vị trí ban đầu của nó sau 10 chu kì thời gian trôi qua (lúc t = 10) cho thấy nó đã lan truyền với tốc độ đúng (c = 1).

wave1da

Hình. Sự lan truyền của xung sóng Gauss 1 chiều. Tính toán số trị được thực hiện với c = 1, δt = 5 × 10−3, và N = 100. Đường liền nét biểu thị Ex lúc t = 0, đường gạch đứt biểu thị Ex lúc t = 10, còn đường chấm chấm (bị đường đứt che khuất) biểu thị By lúc t = 10. [fwave1da]

Hình [fwave1db] cho thấy một ví dụ tính toán thứ ba trong đó dùng đoạn chương trình trên để truyền một xung dạng sóng vuông. Lưu ý rằng lần này chương trình thực hiện rất kém, vì các dao động phi lý đã phát sinh ở những cạnh bên dựng đứng của hình sóng. Như đã thảo luận ở Mục [Sai phân ngược dòng], những dao động như vậy có thể được triệt tiêu bằng cách dùng một lược đồ sai phân ngược dòng, vốn trong trường hợp này nghĩa là sai phân không gian phải bị kéo lệnh về phía hướng truyền sóng. Không may là các lược đồ ngược dòng dạng hiện đơn giản đều chịu ràng buộc CFL,

δt < δz / c, 

và cũng có xu hướng phân tán cao độ.

waveb

Hình. Sự lan truyền của một xung dạng sóng vuông 1 chiều. Tính toán số trị được thực hiện với c = 1, δt = 5 × 10−3, và N = 100. Đường chấm chấm biểu thị Ex lúc t = 0. 0, còn đường liền nét biểu thị Ex lúc t = 0. 1.[fwave1db]

Khoang cộng hưởng hai chiều

Hình [frescav] cho thấy một khoang cộng hưởng 2 chiều gồm một kênh hình hộp chữ nhật có kích thước Lx × Ly. Giả sử rằng các vách kênh chạy dọc theo hai trục xy. Ta sẽ khuấy động kênh này theo cách nhân tạo bằng việc tác động theo phương z một dòng điện xoay chiều có tần số f, vốn có cùng cấu trúc không gian như mode mà ta quan tâm. Ta hãy tính các trường điện từ được khuấy động trong kênh bởi kiểu dòng điện trên.

rescav

Hình. Khoang cộng hưởng hai chiều.[frescav]

Điện trường và từ trường trong kênh có thể lần lượt được viết là E = [0, 0, Ez(x, y, t)], và B = [Bx(x, y, t), By(x, y, t), 0]. Theo phương trình Maxwell ta có

\begin{array}{ccc}  \partial H_x/\partial t + c\,\partial E_z / \partial y &=& 0,\\  \partial H_y/\partial t + c\,\partial E_z / \partial x &=& 0,\\  \partial E_z/\partial t + c\,\partial H_y / \partial x  + c\,\partial H_x / \partial y &=& J_z,\\  \end{array}           [rc1] & [rc3]

trong đó c là vận tốc ánh sáng, Hx = cBx, Hy =  − cBy, và Jz =  − μ0c2jz. Lưu ý rằng hệ phương trình trên có dạng ba phương trình chuyển tải có liên hệ nhau với một số hạng nguồn. Các điều kiện biên gồm có điện trường theo phương tiếp tuyến và từ trường theo phương pháp tuyến đều phải bằng không tại vách dẫn điện. Từ đó suy ra

Ez = 0                       [rcbc1]

ở tất cả các vách (đặt tại x = 0, Lxy = 0, Ly),

Hx = ∂Hy / ∂x = 0

tại x = 0, Lx, và

Hy = ∂Hx / ∂y = 0                   [rcbc2]

tại y = 0, Ly. Sau cùng, kiểu dòng điện được chuẩn hóa ứng với mode (m, n) có dạng

Jz(x, y, t) = J0 sin(mπx / Lx) sin(nπy / Ly) sin(2 πft). 

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, ⋯. Mặt khác, theo phương x, ta làm rời rạc trên lưới đều xi = iδx, với i = 0, I, trong đó δx = Lx / I. Sau cùng, theo phương y, ta làm rời rạc trên lưới đều yj = jδy, với j = 0, J, trong đó δy = Ly / J. Dùng một lược đồ sai phân Crank-Nicolson theo thời gian giống như đã đề cập ở Mục [Lược đồ Crank-Nicolson] và [PT sóng 1 chiều], các PT ([rc1])–([rc3]) cho ta

\begin{array}{rcl}  \displaystyle{\frac{(H_x)_{i,j}^{n+1} - (H_x)_{i,j}^{n}}{\delta t}  +\frac{c}{2}\left( \frac{\partial E_z} {\partial y} \right)^{n+1}_{i,j}  +\frac{c}{2}\left( \frac{\partial E_z} {\partial y} \right)^{n}_{i,j}} &=& 0,\\  \displaystyle{\frac{(H_y)_{i,j}^{n+1} - (H_y)_{i,j}^{n}}{\delta t}  +\frac{c}{2}\left( \frac{\partial E_z} {\partial x} \right)^{n+1}_{i,j}  +\frac{c}{2}\left( \frac{\partial E_z} {\partial x} \right)^{n}_{i,j}} &=& 0,\\  \displaystyle{\frac{(E_z)_{i,j}^{n+1} - (E_z)_{i,j}^{n}}{\delta t}  +\frac{c}{2}\left( \frac{\partial H_y} {\partial x} \right)^{n+1}_{i,j}  +\frac{c}{2}\left( \frac{\partial H_y} {\partial x} \right)^{n}_{i,j}} && \\  \displaystyle{+\frac{c}{2}\left( \frac{\partial H_x} {\partial y} \right)^{n+1}_{i,j}+\frac{c}{2}\left( \frac{\partial H_x} {\partial y} \right)^{n}_{i,j}} &=& (J_z)_{i,j}^n,\\  \end{array}                [xrc1] & [xrc3]

trong đó (Hx)i, jn ≡ Hx(xi, yj, tn), v.v.

Dùng phương pháp Fourier, ta viết

\begin{array}{ccc}  (E_z)_{i,j}^n &=& \sum_{i'=0,I}^{j'=0,J}\hat{E}_{i',j'}^n\,\sin(i\,i'\,\pi/I)\,\sin(j\,j'\,\pi/J),\\  (H_x)_{i,j}^n &=& \sum_{i'=0,I}^{j'=0,J} \hat{X}_{i',j'}^n\,\sin(i\,i'\,\pi/I)\,\cos(j\,j'\,\pi/J),\\  (H_y)_{i,j}^n &=& \sum_{i'=0,I}^{j'=0,J} \hat{Y}_{i',j'}^n\,\cos(i\,i'\,\pi/I)\,\sin(j\,j'\,\pi/J),\\  (J_z)_{i,j}^n &=& \sum_{i'=0,I}^{j'=0,J} \hat{J}_{i',j'}^n\,\sin(i\,i'\,\pi/I)\,\sin(j\,j'\,\pi/J). \\  \end{array}

Những phương trình trên tự thỏa mãn các điều kiện biên ([rcbc1])–([rcbc2]). Các PT ([xrc1])–([xrc3]) cho ta

\begin{array}{ccc}  \hat{X}_{i,j}^{n+1} - \hat{X}_{i,j}^n + j\,D_y\,(\hat{E}_{i,j}^{n+1}+\hat{E}_{i,j}^{n})&=&0,\\  \hat{Y}_{i,j}^{n+1} - \hat{Y}_{i,j}^n + i\,D_x\,(\hat{E}_{i,j}^{n+1}+\hat{E}_{i,j}^{n})&=&0,\\  \hat{E}_{i,j}^{n+1} - \hat{E}_{i,j}^n - i\,D_x\,(\hat{Y}_{i,j}^{n+1}+\hat{Y}_{i,j}^{n}) - j\,D_j\,(\hat{X}_{i,j}^{n+1}+\hat{X}_{i,j}^{n})&=&\delta t\,\hat{J}^{n}_{i,j},\\  \end{array}

với i = 0, Ij = 0, J, trong đó Dx = πcδt / (2 Lx)Dy = πcδt / (2 Ly). Từ đó suy ra

\begin{array}{ccl}  \hat{E}^{n+1}_{i,j} &=& \displaystyle{ \frac{(1-i^2\,D_x^2-j^2\,D_y^2)\,\hat{E}_{i,j}^n + 2\,j\,D_y\, \hat{X}_{i,j}^n+ 2\,i\,D_x \,\hat{Y}_{i,j}^n +\delta t \,\hat{J}^n_{i,j}} {1 + i^2\,D_x^2+j^2\,D_y^2} } \\  \hat{X}_{i,j}^{n+1}&=& \hat{X}_{i,j}^n - j\,D_y\,(\hat{E}_{i,j}^{n+1}+\hat{E}_{i,j}^{n}),\\  \hat{Y}_{i,j}^{n+1}&=& \hat{Y}_{i,j}^n - i\,D_x\,(\hat{E}_{i,j}^{n+1}+\hat{E}_{i,j}^{n}).  \end{array}                           [xstepa] & [xstepb]

Đoạn chương trình dưới đây giải phương trình sóng 2 chiều trong kênh cộng hưởng bằng lược đồ Crank-Nicolson nêu trên. Đầu tiên, chương trình biến đổi Fourier các đại lượng Hx, Hy, Ez, và Jz theo cả 2 phương xy, sau đó chọn bước thời gian theo các PT ([xstepa])–([xstepb]), rồi tái lập Hx, Hy, và Ez qua biến đổi Fourier kép ngược.

// Wave2D.cpp 

// Function to evolve 2-d wave equation: 

//  d H_x / dt + c d E_z / dy = 0 

//  d H_y / dt + c d E_z / dx = 0 

//  d E_z / dt + c d H_y / dx + c d H_x / dy = J_z 

// in region 0 < x < L_x and 0 < y < L_y 

// Boundary conditions: 

//  E_z(0, y) = E_z(L_x, y) = E_z(x, 0) = E_z(x, L_y) = 0 

//  H_x(0, y) = H_x(L_x, y) = d H_y(0, y) / dx =  d H_y(L_x, y) / dx = 0 

//  H_y(x, 0) = H_y(x, L_y) = d H_x(x, 0) / dy =  d H_x(x, L_y) / dy = 0 

// Matrices Hx, Hy, Ez, Jz assumed to be of extent I+1, J+1. 
// Now, (i,j)th elements of matrices correspond to 

//  x_i = i * dx    i=0,I 

//  y_j = j * dy    j=0,J 

// Here, dx = L_x / I is grid spacing in x-direction, 

// and dy = L_y / J is grid spacing in x-direction. 

// Now, Dx = pi c dt / (2 L_x) and Dy = pi c dt / (2 L_y), 

// where dt is time-step. 

// Uses Crank-Nicolson scheme. 

#include <blitz/array.h> 

using namespace blitz; 

void fft_forward_cos (Array<double,1> f, Array<double,1>& F); 
void fft_backward_cos (Array<double,1> F, Array<double,1>& f); 
void fft_forward_sin (Array<double,1> f, Array<double,1>& F); 
void fft_backward_sin (Array<double,1> F, Array<double,1>& f); 

void Wave2D (Array<double,2>& Hx, Array<double,2>& Hy, Array<double,2>& Ez, 
             Array<double,2> Jz, double Dx, double Dy, double dt) 
{ 
  // Find I and J. Declare local arrays 
  int I = Hx.extent(0) - 1; 
  int J = Hx.extent(1) - 1; 
  Array<double,2> X(I+1, J+1), XX(I+1, J+1), XXX(I+1, J+1); 
  Array<double,2> Y(I+1, J+1), YY(I+1, J+1), YYY(I+1, J+1);  
  Array<double,2> E(I+1, J+1), EE(I+1, J+1), EEE(I+1, J+1); 
  Array<double,2> K(I+1, J+1), KK(I+1, J+1); 

  // Fourier transform solution in x-direction 
  for (int j = 0; j <= J; j++) 
    { 
      Array<double,1> In(I+1), Out(I+1); 

      // Fourier transform Hx 
      for (int i = 0; i <= I; i++) In(i) = Hx(i, j); 
      fft_forward_sin (In, Out); 
      for (int i = 0; i <= I; i++) X(i, j) = Out(i); 

      // Fourier transform Hy 
      for (int i = 0; i <= I; i++) In(i) = Hy(i, j); 
      fft_forward_cos (In, Out); 
      for (int i = 0; i <= I; i++) Y(i, j) = Out(i); 

      // Fourier transform Ez 
      for (int i = 0; i <= I; i++) In(i) = Ez(i, j); 
      fft_forward_sin (In, Out); 
      for (int i = 0; i <= I; i++) E(i, j) = Out(i); 

      // Fourier transform Jz 
      for (int i = 0; i <= I; i++) In(i) = dt * Jz(i, j); 
      fft_forward_sin (In, Out); 
      for (int i = 0; i <= I; i++) K(i, j) = Out(i); 
    } 

  // Fourier transform solution in y-direction 
  for (int i = 0; i <= I; i++) 
    { 
      Array<double,1> In(J+1), Out(J+1); 

      // Fourier transform Hx 
      for (int j = 0; j <= J; j++) In(j) = X(i, j); 
      fft_forward_cos (In, Out); 
      for (int j = 0; j <= J; j++) XX(i, j) = Out(j); 

      // Fourier transform Hy 
      for (int j = 0; j <= J; j++) In(j) = Y(i, j); 
      fft_forward_sin (In, Out); 
      for (int j = 0; j <= J; j++) YY(i, j) = Out(j); 

      // Fourier transform Ez 
      for (int j = 0; j <= J; j++) In(j) = E(i, j); 
      fft_forward_sin (In, Out); 
      for (int j = 0; j <= J; j++) EE(i, j) = Out(j); 

      // Fourier transform Jz 
      for (int j = 0; j <= J; j++) In(j) = K(i, j); 
      fft_forward_sin (In, Out); 
      for (int j = 0; j <= J; j++) KK(i, j) = Out(j); 
    } 

  // Evolve XX, YY, and EE 
  for (int i = 0; i <= I; i++) 
    for (int j = 0; j <= J; j++) 
      { 
        double x = double (i) * Dx; 
        double y = double (j) * Dy; 
        double fp = 1. + x*x + y*y; 
        double fm = 1. - x*x - y*y; 

        EEE(i, j) = fm * EE(i, j) + 2. * y * XX(i, j) + 
          2. * x * YY(i,j) + KK(i, j); 
        EEE(i, j) /= fp; 
        XXX(i, j) = XX(i, j) - y * (EEE(i, j) + EE(i, j)); 
        YYY(i, j) = YY(i, j) - x * (EEE(i, j) + EE(i, j));   
      } 

  // Reconstruct solution via inverse Fourier transform in y-direction 
  for (int i = 0; i <= I; i++) 
    { 
      Array<double,1> In(J+1), Out(J+1); 

      // Reconstruct Hx 
      for (int j = 0; j <= J; j++) In(j) = XXX(i, j); 
      fft_backward_cos (In, Out); 
      for (int j = 0; j <= J; j++) X(i, j) = Out(j); 

      // Reconstruct Hy 
      for (int j = 0; j <= J; j++) In(j) = YYY(i, j); 
      fft_backward_sin (In, Out); 
      for (int j = 0; j <= J; j++) Y(i, j) = Out(j); 

      // Reconstruct Ez 
      for (int j = 0; j <= J; j++) In(j) = EEE(i, j); 
      fft_backward_sin (In, Out); 
      for (int j = 0; j <= J; j++) E(i, j) = Out(j); 
    }  

  // Reconstruct solution via inverse Fourier transform in x-direction 
  for (int j = 0; j <= J; j++) 
    { 
      Array<double,1> In(I+1), Out(I+1); 

      // Reconstruct Hx 
      for (int i = 0; i <= I; i++) In(i) = X(i, j); 
      fft_backward_sin (In, Out); 
      for (int i = 0; i <= I; i++) Hx(i, j) = Out(i); 

      // Reconstruct Hy 
      for (int i = 0; i <= I; i++) In(i) = Y(i, j); 
      fft_backward_cos (In, Out); 
      for (int i = 0; i <= I; i++) Hy(i, j) = Out(i); 

      // Reconstruct Ez 
      for (int i = 0; i <= I; i++) In(i) = E(i, j); 
      fft_backward_sin (In, Out); 
      for (int i = 0; i <= I; i++) Ez(i, j) = Out(i); 
    } 
}

Việc tính toán số trị như dưới đây được thực hiện bằng đoạn chương trình trên. Các trường điện từ Hx, Hy, và Ez đều được khởi tạo bằng 0 tại mọi điểm lúc t = 0. Hình [fcavity1] biểu diễn biên độ cực đại của Ez theo tần số, f, của một phân bố dòng điện m = 1 / n = 1. Có thể thấy rằng rõ ràng có cộng hưởng tại f ≃ 0. 7.

cavity1

Hình. Các sóng điện từ trong một kênh cộng hưởng 2 chiều. Biên độ cực đại của Ez(x = 0. 5, y = 0. 5) từ t = 0 đến t = 20 theo tần số nguồn, f. Tính toán số trị được thực hiện với m = 1, n = 1, Lx = 1, Ly = 1, c = 1, I = J = 32, và δt = 10−2.[fcavity1]

Các Hình [fcavity2] và [fcavity3] minh họa sự biến đổi theo thời gian điển hình của Ez, Hx, và Hy lần lượt trong trường hợp không cộng hưởng và cộng hưởng. Với trường hợp không cộng hưởng, các vết có dạng mẫu can thiệp giữa phản hồi trực tiếp, vốn dao động với tần số nguồn f, và phản hồi tức thời, vốn dao động với tần số tự nhiên f0 của kênh. Lưu ý rằng các dao động tức thời không bao giờ tắt, vì trong bài toán này ta không xét sự tiêu tán năng lượng. Cũng cần nói thêm là ta dễ thấy được

f_0 = \frac{c}{2} \sqrt{ \displaystyle{ \frac{1}{(n\,L_x)^2}+\frac{1}{(m\,L_y)^2} } }.

Điều này dẫn đến f_0=1/\sqrt{2} = 0.7071 với n = m = Lx = Ly = c = 1, vốn tương ứng hoàn toàn với tần số cộng hưởng ở Hình [fcavity1]. Trong trường hợp cộng hưởng, các vết có dạng sóng với biên độ ngày càng tăng dần và dao động với tần số tự nhiên f0.

cavity2

Hình. Các sóng điện từ trong kênh cộng hưởng 2 chiều. Các vết theo thời gian của Ez(x = 0. 5, y = 0. 5) (nét liền), Hx(x = 0. 5, y = 0. 0) (nét đứt), và Hy(x = 0. 0, y = 0. 5) (nét chấm—bị che khuất bởi nét đứt). Tính toán số trị được thực hiện với m = 1, n = 1, Lx = 1, Ly = 1, c = 1, I = J = 32, δt = 10−2, và f = 0. 6.[fcavity2]

cavity3

Hình. Các sóng điện từ trong kênh cộng hưởng 2 chiều. Các vết theo thời gian của Ez(x = 0. 5, y = 0. 5) (nét liền), Hx(x = 0. 5, y = 0. 0) (nét đứt), và Hy(x = 0. 0, y = 0. 5) (nét chấm—bị che khuất bởi nét đứt). Tính toán số trị được thực hiện với m = 1, n = 1, Lx = 1, Ly = 1, c = 1, I = J = 32, δt = 10−2, và f = 0. 7071.[fcavity3]

Sau cùng, các Hình [fcavity4] và [fcavity5] minh họa sự thay đổi trong không gian của các trường điện từ được gây ra trong kênh khi m = 1n = 1.

cavity4

Hình. Các sóng điện từ trong kênh cộng hưởng 2 chiều. Sự thay đổi trong không gian của Ez (nét liền) và Hy (nét đứt) theo phương x tại t = 20y = 0. 5. Tính toán số trị được thực hiện với m = 1, n = 1, Lx = 1, Ly = 1, c = 1, I = J = 32, δt = 10 − 2, và f = 0. 7071.[fcavity4]

cavity5

Hình. Các sóng điện từ trong kênh cộng hưởng 2 chiều. Sự thay đổi trong không gian của Ez (nét liền) và Hx (nét đứt) theo phương y tại t = 20x = 0. 5. Tính toán số trị được thực hiện với m = 1, n = 1, Lx = 1, Ly = 1, c = 1, I = J = 32, δt = 10−2, và f = 0. 7071.[fcavity5]

Advertisements

%(count) bình luận

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

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

  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 Đă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