AEP 4830 HW6 Boundary Value Problem and Relaxation

Finite Difference Algorithm

We have used Runge-Kutta method to solve ODEs with initial conditions; however, many of the physical systems cannot be described by ODEs. The partial differential equations with boundary conditions are inevitable for physical systems such as time-independent Schrodinger equations and electromagnetic fields. Those equations involve partial differential in spatial coordinates with specified boundaries and/or sources. When trying to solve such system numerically, the finite difference algorithm is a straightforward and easy-to-program method.
We will use the 2D Laplace equation as an example to demonstrate the algorithm.

\begin{equation} \nabla^{2}\phi(x,y)=0\\ \end{equation}

The finite difference algorithm requires to mesh-grid the space of interest and the value on each grid point \((i,j)\) represents the \(\phi(x=ih_{1},y=jh_{2})=\phi_{ij}\) where \(h_{1}\) and \(h_{2}\) are the spacing in \(x\) and \(y\) axes. It’s always required to set up another ”flag” grid which determines whether certain combination of \((i,j)\) is on the boundary or source (electrodes). Next, based on the finite difference of adjacent values to \(\phi_{ij}\) and spacing, the second partial differential can be approximated. We can solve for \(\phi_{ij}\) in terms of its neighbors for the first iteration.(W. H. Press 2007)

\begin{equation} \phi_{ij}^{FD}=\frac{1}{4}\left(\phi_{i-1,j}+\phi_{i+1,j}+\phi_{i,j+1}+\phi_{i,j-1}\right)\\ \end{equation}

The flag grid is like a mask, leaving all the specified boundary values unchanged during iterations. Similar to Runge-Kutta method, we evaluate the maximum error in each iteration.

\begin{equation} \Delta\phi_{ij}=max|\phi_{ij}^{FD}-\phi_{i,j}|\\ \end{equation}

If the maximum error is larger than the specified tolerance \(\Delta\phi_{ij}>\epsilon_{max}\), we should pass the new values of \(\phi_{ij}\) to Eq. (2) and repeat the calculation until the error is below the tolerance or the number of iterations exceeds certain upper bound. The new \(\phi_{ij}^{new}\) can be the following form, which is the weighed average of \(\phi_{ij}\) and \(\phi_{ij}^{FD}\).

\begin{equation} \phi_{ij}^{new}=(1-\omega)\phi_{ij}+\omega\phi_{ij}^{FD}\\ \end{equation}

where \(0<\omega<2\). If \(1<\omega<2\), the algorithm converges more quickly and this is known as the successive over-relaxation algorithm. The regime where \(0<\omega<1\) corresponds to the under-relaxation where the solution converges relatively slowly.
Here, we will solve the Laplace equation in a cylinder with three disks of inner and outer radius. Instead of utilizing Eq. (2), we will use the finite difference in cylindrical coordinates and assume azimuthal symmetry.

\begin{equation} \phi_{ij}^{FD}=\frac{1}{4}\left(\phi_{i-1,j}+\phi_{i+1,j}+\phi_{i,j+1}+\phi_{i,j-1}\right)+\frac{1}{8j}\left(\phi_{i,j+1}-\phi_{i,j-1}\right)\\ \end{equation}

,where \(z=ih\) and \(r=jh\). The above equation applies to points where \(r\neq 0\). On the central axis where \(r=0\), we should have

\begin{equation} \phi_{ij}^{FD}=\frac{1}{6}\left(4\phi_{i,j=1}+\phi_{i+1,j=0}+\phi_{i-1,j=0}\right)\\ \end{equation}

In order to use finite difference algorithm, we need to set up matrices in c++. We used arrayt.h to generate and do the calculations for matrices.


First, we investigated the powerful role of \(\omega\) here. We made the grid spacing \(h=0.1\) and error tolerance \(\epsilon_{max}=2\) and calculated the number of iterations for the error to go below the tolerance. The results are shown in the following table.

\(\omega\) 1.0 1.2 1.4 1.6 1.8 1.9
Number of iterations 655 511 431 371 359 1077

The best choice of \(\omega\) seems to be around \(\omega=1.8\) and we would use this over-relaxation value for the following calculation.


For solving the electric potential inside the cylinder with three disk electrodes and boundaries, the grid spacing is the same \(h=0.1\) while the error tolerance is much smaller, \(\epsilon_{max}=0.01\). This one percent error would increase the number of iterations significantly to \(N=4231\). The contour plot of the numerical solution in space is shown in Fig. 1, where the red line represents the middle disk electrode held at \(\psi=2000V\) and the black line and rectangle are the two grounded electrodes \(\psi=0V\). The electric potential contour is distorted by the two grounded electrodes, compressing the potential lines in the z direction. The potential at about \(r=10mm\) is symmetric while that closer to the central axis is slightly asymmetric due to the mismatch in inner radii of the two grounded electrodes. The electric potential is mainly restricted among the disks and near the origin. Outside of those disks, the potential is negligibly zero everywhere. The configuration is seen in the electron guns of the electron microscopes to focus electrons by adjusting the fixed voltage of the center electrode.

The electric potential in space where the value of each contour line is labeled. The middle electrode is held at \(\phi=2000V\) while the other black two are grounded. The boundaries are all grounded. Notice the asymmetry if the potential lines due to the gap difference and the inner radii mismatch near the central axis.

Since we now have solved the potential everywhere in space, we may as well calculate the electric fields. The elect