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 \begin{equation} \omega_{nm}=\sqrt{-k^2}=\frac{\xi_n^{(m)}}{R}\;. \end{equation}
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.
