2-d problem with Dirichlet boundary conditions

Let us consider the solution of Poisson's equation in two dimensions. Suppose that

$$\frac{\partial^2 u(x,y)}{\partial x^2}+\frac{\partial^2 u(x,y)}{\partial y^2} = v(x,y),$$ (146)

for $x_l\leq x\leq x_h$, and $0\leq y \leq L$. By direct analogy with our previous method of solution in the 1-d case, we could discretize the above 2-d problem using a second-order, central difference scheme in both the $x$- and $y$-directions. Unfortunately, such a discretization scheme yields a set of equations which cannot be reduced to a simple tridiagonal matrix equation. In fact, all of the efficient numerical algorithms for solving this type of problem are iterative in nature. For instance, the Jacobi method, the Gauss-Seidel method, the successive over-relaxation method, and the multi-grid method.³⁴Regrettably, unless such iteration methods are extremely sophisticated (e.g., the multi-grid method), and, hence, beyond the scope of this course, they tend to converge very poorly. In the following, rather than discuss iterative methods which do not work very well, we shall instead discuss a non-iterative method which works effectively for a restricted set of problems. The method in question is termed a spectral method, since it involves expanding $u$ and $v$ as truncated Fourier series in the $y$-direction.

Suppose that $u(x,y)$ satisfies mixed boundary conditions in the $x$-direction: i.e.,

$$\alpha_l\,u(x,y) + \beta_l\,\frac{\partial u(x,y)}{\partial x} = \gamma_l(y),$$ (147)

at $x=x_l$, and

$$\alpha_h \,u(x,y)+ \beta_h\,\frac{\partial u(x,y)}{\partial x} = \gamma_h(y),$$ (148)

at $x=x_h$. Here, $\alpha_l$, $\beta_l$, etc., are known constants, whereas $\gamma_l$, $\gamma_h$ are known functions of $y$. Furthermore, suppose that $u(x,y)$ satisfies the following simple Dirichlet boundary conditions in the $y$-direction:

$$u(x,0) = u(x,L) = 0.$$ (149)

Note that, since $u(x,y)$ is a potential, and, hence, probably undetermined to an arbitrary additive constant, the above boundary conditions are equivalent to demanding that $u$ take the same constant value on both the upper and lower boundaries in the $y$-direction.

Let us write $u(x,y)$ as a Fourier series in the $y$-direction:

$$u(x,y) = \sum_{j=0}^{\infty} U_j(x)\,\sin(j\,\pi\,y/L).$$ (150)

Note that the above expression for $u$ automatically satisfies the boundary conditions in the $y$-direction. The $\sin(j\,\pi\,y/L)$ functions are orthogonal, and form a complete set, in the interval $y=0,L$. In fact,

$$\frac{2}{L}\int_0^L \sin(j\,\pi\,y/L)\,\sin(k\,\pi\,y/L)\,dy = \delta_{jk}.$$ (151)

Thus, we can write the source term as

$$v(x,y) = \sum_{j=0}^{\infty} V_j(x)\,\sin(j\,\pi\,y/L),$$ (152)

where

$$V_j(x) = \frac{2}{L} \int_0^L v(x,y)\,\sin(j\,\pi\,y/L)\,dy.$$ (153)

Furthermore, the boundary conditions in the $x$-direction become

$$\alpha_l\,U_j(x) + \beta_l\,\frac{dU_j(x)}{dx} = \Gamma_{l\,j},$$ (154)

at $x=x_l$, and

$$\alpha_h\,U_j(x)+ \beta_h\,\frac{dU_j(x)}{dx} = \Gamma_{h\,j},$$ (155)

at $x=x_h$, where

$$\Gamma_{l\,j} = \frac{2}{L} \int_0^L \gamma_l(y)\,\sin(j\,\pi\,y/L)\,dy,$$ (156)

etc.

Substituting Eqs. (150) and (152) into Eq. (146), and equating the coefficients of the $\sin(j\,\pi\,y/L)$ (since these functions are orthogonal), we obtain

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

for $j=0,\infty$. Now, we can discretize the problem in the $y$-direction by truncating our Fourier expansion: i.e., by only solving the above equations for $j=0,J$, rather than $j=0,\infty$. This is essentially equivalent to discretization in the $y$-direction on the equally-spaced grid-points $y_j=j\,L/J$. The problem is discretized in the $x$-direction by dividing the domain into equal segments, according to Eq. (114), and approximating $d^2/dx^2$ via the second-order, central difference scheme specified in Eq. (115). Thus, we obtain

$$U_{i-1,j} - (2+j^2\,\kappa^2)\,U_{i,j} + U_{i+1,j} = V_{i,j}\,(\delta x)^2,$$ (158)

for $i=1,N$ and $j=0,J$. Here, $U_{i,j} \equiv U_j(x_i)$, $V_{i,j} \equiv V_j(x_i)$, and $\kappa = \pi\,\delta x/L$. The boundary conditions (154) and (155) discretize to give:

$$U_{0,j} = \frac{\Gamma_{l\,j}\,\delta x-\beta_l\,U_{1,j}}{\alpha_l\,\delta x -\beta_l},$$ (159)

$$U_{N+1,j} = \frac{\Gamma_{h\,j}\,\delta x + \beta_h\,U_{N,j}}{\alpha_h\,\delta x +\beta_h},$$ (160)

for $j=0,J$. Eqs. (158), (159), and (160) constitute a set of $J+1$ uncoupled tridiagonal matrix equations (with one equation for each separate $j$ value). These equations can be inverted, using the algorithm discussed in Sect. 5.4, to give the $U_{i,j}$. Finally, the $u(x_i, y_j)$ values can be reconstructed from Eq. (150). Hence, we have solved the problem.