Accuracy of numerical methods in finding eigenvalues

Question

I have a problem with assesing the accuracy of my numerical calculation. I have a 2nd order ODE. It is an eigenvalue problem of the form: $$ y'' + ay' + \lambda^2y = 0, $$ and the boundary condiations are: $$ y(0) = y(1) = 0. $$

This equation describes a vibrating string, clamped at x=0 and x=1, with a certain mass distribution. I want to be able to calculate the eigenvalues of this equation numerically. I do this by using the Runge-Kutta method to find solutions to the equation with initial conditions: $$ y(0) = 0, y'(0) = 1 $$ with different values of $\lambda$ and then looking for the ones that are zero at x=1.

I terminate my search for eigenvalues when I find a function that has $$y(1) < \delta,$$ where $\delta$ is a predefined constant.

Now the problem is that I'd like to know how accurate my calculations are, i.e. how do I choose the values of $\delta$ and the stepsize used in computing the values of the potential eigenfunction, so that the error in the eigenvalue is less than, say, $ 10^{-3}$?

Answer 1

The exact solution as well as the numerical solution for fixed stepsize $h$ depends analytically on $λ$, at simple eigenvalues there is locally a near linear relationship between $λ$ and the right boundary value.

Thus for a method of order $p$, the numerical eigenvalue for a given stepsize $h$ will lead to a right boundary value of size $O(h^p)$ in the exact solution and thus a distance of the same magnitude to the eigenvalue of the exact solution.

To get an error of magnitude $10^{-3}$ with classical RK4 as integrator one can try $h=0.1$ to $h=0.2$ as step size to get an error of $10^{-4}$ to $10^{-3}$ in the numerical integration and hopefully with a factor close to $1$ the same error magnitude in the eigenvalue. Exact a-priory estimates are difficult, it is easier to estimate the error a-posteriori, for instance using Richardson extrapolation.

---

Note that the higher oszillation for larger $λ$ leads to consider the requirements of the sampling theorem, i.e., the sampling density has to be high enough. Or in other words, for large $λ$ the problem becomes stiff, the error term could be modified to $O(e^{|λ|}h^p)$ to reflect that.

---

This should be no problem for the lowest eigenvalue, thus find the solution with error $10^{-4}$ in the root finder for $h=0.1$ and $h=0.2$ and use the extrapolation formula to estimate the error $$ error = \frac{|λ_{0.2}-λ_{0.1}|}{15} $$ for $λ_{0.1}$.

Answer 2

In principle, you could try to come up with an analytic expression for how the error in step size or $\delta$ affects the answer. I would not bother. You will probably get it wrong. The far, far better method is to vary those numbers and see how much your answer changes. For example, you could change $\delta$ from 0.01 to 0.001 and see how much the answer changes. Similarly, you can vary the step size by factors of 2 or 10.