The 1-d wave equation

Consider a plane polarized electromagnetic wave propagating in vacuo along the $z$-axis. Suppose that the electric and magnetic fields take the form ${\bf E} = [E_x(z,t),0,0]$, and ${\bf B} = [0,B_y(z,t),0]$. Now, according to Maxwell's equations,

$$\frac{\partial E_x}{\partial t} + c\,\frac{\partial H_y}{\partial z} = 0,$$ (253)

$$\frac{\partial H_y}{\partial t} + c\,\frac{\partial E_x}{\partial z} = 0,$$ (254)

where $H_y = c\,B_y$, and $c$ is the velocity of light. Note that the above equations take the form of two coupled advection equations. Let us find the numerical solution of these equations in some region $0\leq z \leq L$ which is bounded by perfectly conducting walls at $z=0$ and $z=L$. Now, both the tangential electric field and the normal magnetic field must be zero at a perfect conductor. It follows that

$$E_x(0,t) = E_x(L,t) = 0.$$ (255)

Moreover, Eq. (253) yields

$$\frac{\partial H_y(0,t)}{\partial z} = \frac{\partial H_y(L,t)}{\partial z} = 0.$$ (256)

We expect Eqs. (253) and (254) to support wave-like solutions which bounce backwards and forwards between the conducting walls.

As before, we discretize in time on the uniform grid $t_n=t_0+n\,\delta t$, for $n=0,1,2,\cdots$. Furthermore, in the $z$-direction, we discretize on the uniform grid $z_i = i\,\delta z$, for $i=0,I$, where $\delta z = L/I$. Adopting a Crank-Nicholson temporal differencing scheme similar to that discussed in Sect. 7.4, Eqs. (253)-(254) yield

$$\frac{(E_x)_i^{n+1}-(E_x)_i^n}{\delta t} +\frac{c}{2}\left(\frac{... ...z}\right)_i^{n+1} + \frac{c}{2}\left(\frac{\partial H_y}{\partial z}\right)_i^n = 0,$$ (257)

$$\frac{(H_y)_i^{n+1}-(H_y)_i^n}{\delta t} +\frac{c}{2}\left(\frac{... ...z}\right)_i^{n+1} + \frac{c}{2}\left(\frac{\partial E_x}{\partial z}\right)_i^n = 0,$$ (258)

where $(E_x)_i^n\equiv E_x(z_i,t_n)$, etc.

Adopting a Fourier approach, we write

$$(E_x)_i^n = \sum_{j=0,I} \hat{E}_j^n\,\sin(i\,j\,\pi/I),$$ (259)

$$(H_y)_i^n = \sum_{j=0,I} \hat{H}_j^n\,\cos(i\,j\,\pi/I),$$ (260)

which automatically satisfies the boundary conditions (255) and (256). Equations (257) and (258) yield

$$\hat{E}_i^{n+1}-\hat{E}_i^n - i\,D\,(\hat{H}_i^{n+1}+\hat{H}_i^n) = 0,$$ (261)

$$\hat{H}_i^{n+1}-\hat{H}_i^n + i\,D\,(\hat{E}_i^{n+1}+\hat{E}_i^n) = 0,$$ (262)

where $D=\pi\,c\,\delta t/(2\,L)$. It follows that

$$\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,$$ (263)

$$\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,$$ (264)

for $i=0,I$.

Let us determine the eigenvalues $\lambda$ of the above linear system, assuming that $(\hat{E}_i^{n+1}, \hat{H}_i^{n+1}) = \lambda\,(\hat{E}_i^n, \hat{H}_i^n)$. After a little analysis, we obtain

$$\lambda^2 - 2\,\frac{1-i^2\,D^2}{1+i^2\,D^2}\,\lambda + 1 = 0.$$ (265)

For all values of $i$, the two roots of the above quadratic are complex, but have modulus unity: i.e., $\vert\lambda\vert=1$. This implies that our differencing scheme is both numerically stable (if $\vert\lambda\vert>1$ then the scheme would be unstable, since the associated eigenfunction would be amplified at each time-step) and non-dispersive (if $\vert\lambda\vert<1$ then all Fourier harmonics would eventually decay away, and, hence, so would the solution). Note that the above calculation is equivalent to von Neumann stability analysis.

The routine listed below solves the 1-d wave equation using the Crank-Nicholson scheme discussed above. The routine first Fourier transforms $E_x$ and $H_y$, takes a time-step using Eqs. (263) and (264), and then reconstructs $E_x$ and $H_y$ via an inverse Fourier transform.

// 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 < z < 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-Nicholson 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 Wave1D (Array<double,1>& Ex, Array<double,1>& Hy, double D)
{
  // Find I. Declare local arrays
  int I = Ex.extent(0) - 1;
  Array<double,1> 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 <= 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);
}

Figure 230 shows an example calculation which uses the above routine to propagate a Gaussian wave-pulse. The initial condition is

$$E_x(z,0) = H_y(z,0) = \exp[-100\,(z-0.5)^2],$$ (266)

and the calculation is performed with $L=1$, $c=1$, $\delta t = 5\times 10^{-3}$, and $N=100$. It can be seen that the pulse moves to the right, reflects off the conducting wall at $z=1$, and then moves to the left. Note, that $E_x=+H_y$ when the wave is far from the conducting walls and moving to the right, whereas $E_x=-H_y$ when the wave is far from the conducting walls and moving to the left.

Figure 80: Propagation of a 1-d Gaussian wave-pulse. Numerical calculation performed using $c=1$, $\delta t = 5\times 10^{-3}$, and $N=100$. The solid curves show $E_x$, whereas the dashed curves show $H_y$. The top-left, top-right, middle-left, middle-right, bottom-left, and bottom-right panels show the solution at $t=0.0$, $0.2$, $0.4$, $0.6$, $0.8$, and $1.0$, respectively.

Figure 253 shows a second example calculation performed with the same parameters as the first. In this calculation, the pulse is allowed to reflect off the conducting walls ten times before returning to its initial position. Note that the pulse amplitude and shape remain constant to a very good approximation during this process. The fact that the pulse returns almost exactly to its initial position after ten time periods have elapsed (i.e., at $t=10$) demonstrates that it is propagating at the correct speed (i.e., $c=1$).

Figure 81: Propagation of a 1-d Gaussian wave-pulse. Numerical calculation performed using $c=1$, $\delta t = 5\times 10^{-3}$, and $N=100$. The solid curve shows $E_x$ at $t=0.$, the dashed curve shows $E_x$ at $t=10.$, and the dotted curve (obscured by the dashed curve) shows $B_y$ at $t=10.$

Figure 82 shows a third example calculation which uses the above listed routine to propagate a square wave-pulse. Note that the routine does a very poor job, since spurious oscillations are generated at the sharp leading and trailing edges of the wave-form. As discussed in Sect. 7.5, such oscillations can only be suppressed by adopting an upwind differencing scheme, which in this case means that the spatial differences must be skewed in the direction from which the wave is propagating. Unfortunately, simple explicit upwind schemes are subject to the CFL constraint,

$$\delta t < \frac{\delta z}{c},$$ (267)

and also tend to be highly dispersive.

Figure 82: Propagation of a 1-d square wave-pulse. Numerical calculation performed using $c=1$, $\delta t = 5\times 10^{-3}$, and $N=100$. The dotted curve shows $E_x$ at $t=0.0$, and the solid curve shows $E_x$ at $t=0.1$.