# Laplacian eigenmodes on a semi-circular region with finite-difference method

The computation of eigenmodes of a semi-circular membrane reduces to the following eigenvalue problem

$$\nabla^2u=k^2u\;,$$

where the region of interest is a semi-circle defined by $r\in[0,1]$ and $\varphi\in[0,\pi]$.

It is appropriate to work in cylindrical coordinates, where the Laplacian is written as

$$\nabla^2u=\frac{\partial^2u}{\partial r^2}+\frac{1}{r}\frac{\partial u}{\partial r}+\frac{1}{r^2}\frac{\partial^2 u}{\partial\varphi^2}.$$

Boundary conditions fix the value of $u$ at the boundary of the semi-circle, where $u=0$.

First, we make a discretization of $u$ with $u_{ij}=u(r_i,\varphi_j)$, where $r_i=(i+\frac{1}{2})h_r$ and $\varphi_j=(j+\frac{1}{2})h_\varphi$ $i,j=0\dots N-1$ and $h_r=1/N$, $h_r=\pi/N$. This is a centered mesh.

We then use a finite difference approximation for the Laplacian and obtain

$$\begin{eqnarray} \nabla^2u&\approx& \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h_r^2}+\frac{1}{ih_r}\frac{u_{i+1,j}-u_{i-1,j}}{2h_r}+ \\ &&+\frac{1}{(ih_r)^2}\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_{\varphi}^2}\\ &=&k^2u_{ij} \end{eqnarray}$$

or

$$\begin{eqnarray} &&u_{i+1,j}\left(1+\frac{1}{2i}\right)+u_{i-1,j}\left(1-\frac{1}{2i}\right)\\ &&+\frac{1}{i^2h_{\varphi}^2}\left(u_{i,j+1}+u_{i,j-1}\right)+u_{i,j}\left(-2-\frac{2}{i^2h_{\varphi}^2}-k^2h_r^2\right)=0\;. \end{eqnarray}$$

Because our mesh is centered we have to make the following replacement in the above equation: $i\to i+\frac{1}{2}$. This replacement also helps us get rid of the coordinate singularity for $i=0$.

Boundary conditions at $\varphi=0,\pi$ and $r=0,1$ can be all handled with the same trick, where we set at the boundary

$$u_{i,j-1}=-u_{i,j}$$ $$u_{i,j+1}=-u_{i,j}$$ $$u_{i-1,j}=-u_{i,j}$$ $$u_{i+1,j}=-u_{i,j}.$$

From $u_{ij}$ we then form a vector $\vec v$ and get a classical eigenvalue problem for a matrix $\rm A$, which is carefully formed from above equations $${\rm A}{\vec v}=k^2h_r^2{\vec v}\;.$$

The matrix is an unsymmetric real matrix and eigenvalues and eigenvectors can be obtained with a routine dgeev from LAPACK.

Analytical solutions can easily be obtained by the method of separation of variables

$$u(r,\varphi)=R(r)\Phi(\varphi)\;.$$

They are

$$u(r,\varphi)_{nm}=\sin(n\varphi)J_n\left(\frac{\xi_n^{(m)}}{R}r\right)\;,$$ where $J_n$ is a (cylindrical) Bessel function of the first kind of order $n$ and $\xi_n^{(m)}$ is the $m$-th zero of $J_n$.

Eigenvalues and frequencies are $$\omega_{nm}=\sqrt{-k^2}=\frac{\xi_n^{(m)}}{R}\;.$$

My problem is that the numerical solution obtained by the described procedure does not match the analytical one. The difference is around $r=0$, so this boundary probably isn't correctly accounted for. My results can be seen in the following plots of both analytical and numerical solutions.

Here is the plot of the analytical solution for the first eigenfunction:

The following plot shows the comparison of numerical results for three different discretizations, as far as my computational resources allow me to go.

The next plot shows the $N^2$ dependence of differences between numerical and analytical solutions, normalized as $L^2(\vec{u}-\vec{u}_{analytical})/N^2$ in logarithmic scale. The slope of the linear regression is $0.5$, which means that absolute error is linearly decreasing with $N$. This linear accuracy is not surprising, since the boundary conditions were satisfied to first order only, though other difference approximations were second order.

• I'm a little confused. You're solving the discrete eigenproblem corresponding to $A$, correct? Its eigenvalues and eigenvectors should be different than those of the continuous problem, due to discretization error, so I would expect that the results would not be the same. As $h_{r}, h_{\varphi} \rightarrow 0$, the spectrum of the discrete problem should approach the (infinite) spectrum of the continuous problem, and the eigenvectors of $A$ should approach the values of the eigenfunctions of the continuous problem at the grid points of your simulation. – Geoff Oxberry Jan 23 '12 at 17:12
• That's true. But as you can see from the posted solution (where nb. of points $N=60$ was used) that the numerical solution isn't approaching the analytical one. The derivatives around $r=0$ seem to be too small. I think some boundary condition at $r=0$ should be taken in account. But the condition $u(r=0,\varphi)=0$, if implemented in same way as for other boundaries (as described above) doesn't change anything. – liberias Jan 23 '12 at 18:18
• @liberias you can't say anything about convergence from just one result with one value of $N$. Why not try taking $N$ larger and see if your result improves? Better yet, do a real convergence test. – David Ketcheson Jan 23 '12 at 20:50
• @liberias Nice! Can you compute the (appropriately scaled) norm of the difference between the numerical and analytical solutions for the different values of $N$? It would be useful to see how quickly it is decreasing. – David Ketcheson Jan 24 '12 at 16:39
• @David Is this simple convergence test what you meant? Is this (linear) convergence fast enough? After all we used second order approximation for the first derivative. – liberias Jan 28 '12 at 13:16

In the interest of having at least some answer, it seems like the analysis you have above is correct. First-order boundary conditions do affect the accuracy of a second-order discretization (if I'm remembering LeVeque's book on finite difference approximations correctly; DavidKetcheson can help me out with that), so the solutions of the discrete problem appear to be converging to the solution of the analytical problem at the appropriate rate.

• Gilbert Strang also mentions in his lectures that first order boundary conditions propagate throughout the solution. Actually I am wondering if this kind of behaviour near origin is typical for cylindrical coordinates and possibly originates in coordinate singularity. – liberias Feb 9 '12 at 13:50

Having skimmed over most of the analysis and looked at the pretty pictures, they do tell a convincing story: the analytical solution does not have any circularity to it, but your numeric solution does; the central region is curved.

So I would focus squarely on your finite difference scheme. I can't see what's wrong yet, but perhaps some part of the angular dependence of the coordinates is not quite completely taken care of.

I doubt it is the boundary condition; both solutions seem to have zero at all the boundaries. It could conceivably be the nature of the condition though, being a derivative rather than value condition.

I do not understand the boundary conditions.

• Why the minus?
• These four conditions holds for $r=0$ AND $r=1$ AND $\varphi=0$ AND $\varphi=\pi$? Why?

Could you give more details, please?

• Why the minus? Suppose you have a linear function $u(x)$ which has a value $u_1$ at point $x=-a$ and value $-u_1$ at point $x=a$. Then the function's value at point $x=0$ is equal $0$. Why four boundary conditions? When working with coordinates $r$ and $\varphi$ on a square domain there are four boundaries ($\varphi=0,\pi$ and $r=0,1$) on which you need to specify some boundary conditions. – liberias Jul 2 '12 at 9:18