Numerically solving a BVP for a 2nd order ODE with mixed Dirichlet and Neumann boundary conditions with the finite difference method

Question

I would like to solve a simple BVP for a second order differential equation on the domain $I=[0,\pi/2]$ using the finite difference method. $$ u''(x)+u(x)=0\\ u(0)=0\\ u'(\pi/2)=0 $$ By construction, $u(x) = c \sin(x)$ is a solution for any $c$. To use the finite difference method, we construct an $n\times n$ matrix of the form $$ A = \frac{1}{\Delta x^2}\begin{bmatrix} -2&1&0&0&\cdots&0\\ 1&-2&1&0&\cdots&0\\ \ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\ 0&\cdots&0&1&-2&1\\ 0&\cdots&0&0&2&-2 \end{bmatrix}. $$ Here in the first row we are assuming $u_0=0$ for the zero Dirichlet condition, and we are assuming $u_n=u_{n+1}$ so that $u'(x_n)=0$. Then we solve the discretized differential equation, which becomes a linear system, $$ A\vec{u}+\vec{u}=0\\ (A+I)\vec{u}=0, $$ but since $(A+I)$ is invertible, the only solution to this equation is $\vec{u}=0$. How can I recover the non-trivial solutions?

Answer 1

The point is that the number $\pi/2$, that is to say the domain size in which these boundary conditions result in an eigenvalue of $-1$, is part of the solution of the problem. It will come with its own numerical error. This formulation means you are asked to find $L$ such that $u''+u=0,u(0)=0,u'(L)=0$ has a nontrivial solution. This $L$ is $n\Delta x$, so you want $A+I$ to be singular, or equivalently you want $\Delta x^2 A + \Delta x^2 I$ to be singular. So you want $\Delta x$ to be $\sqrt{-\lambda}$ for an eigenvalue $\lambda$ of the matrix $\Delta x^2 A$ (which depends only on $n$ and not on $\Delta x$).

If you don't like the idea of the domain size as a variable, then you can hold the domain size fixed at $\pi/2$ and allow for the numerical eigenvalue to differ from $-1$ instead.

Either way as you observed you need to leave some slack, because if the domain is size $\pi/2$ then the eigenvalue you get is never exactly $-1$, and conversely if you set the eigenvalue to be exactly $-1$ then the domain size you get is not exactly $\pi/2$.