Accuracy of numerical methods in finding eigenvalues

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}$?

• personally, I would use a variation on MatLab's "fzero" that has adaptive step-sizes. As long as the initial value is in the neighborhood, and the initial step-sizes are not too big, the accuracy is arbitrary - which means you can set it to $10^{-10}$ if you like. – EngrStudent Sep 14 '15 at 16:57
• It seems like this will only find solutions satisfying $y'(0)=1$, which isn't part of the problem specification. – David Ketcheson Sep 14 '15 at 17:46
• The bigger problem is that with the procedure you use, you can only find real eigenvalues. But, for $a\neq 0$, your eigenvalues will have a nonzero imaginary component. – Wolfgang Bangerth Sep 14 '15 at 18:41
• In the preceding comnent, the OP write "as a part of the excercise, I cannot use the inbuilt functions". Your question should have begun "I have a homework problem I want you to do for me". – Mark L. Stone Sep 15 '15 at 3:14
• There are a dozen ways to find a zero. Bisection is cheap. Fzero uses a secant-Newton method. Any decent numerical recipes text has chapters dedicated to the topic. – EngrStudent Sep 15 '15 at 3:31

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}$.

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.