6
$\begingroup$

I want to find the roots for $\kappa$ for the equation

$$\sqrt{\alpha - 1} \cos{\left (\frac{\sqrt{2} \sqrt{\alpha - 1}}{2 \sqrt{\epsilon}} \right )} \cosh{\left (\frac{\sqrt{2} \sqrt{\alpha + 1}}{2 \sqrt{\epsilon}} \right )} - \sqrt{\alpha - 1} \\ -\frac{1}{\sqrt{\alpha + 1}} \sin{\left (\frac{\sqrt{2} \sqrt{\alpha - 1}}{2 \sqrt{\epsilon}} \right )} \sinh{\left (\frac{\sqrt{2} \sqrt{\alpha + 1}}{2 \sqrt{\epsilon}} \right )} = 0 \enspace ,$$ where $\alpha=\sqrt{1 + 4\epsilon\kappa^2}$. This equation has infinite roots, but I am interested in the first $N$ of them.

One option to solve this problem is to use Newton method, the problem is to define the initial points, since the values of the function can be quite high, as can be seen below

enter image description here

This problem comes from finding the eigenvalues of

$$\left[\frac{d^2 u}{ds^2} - \epsilon \frac{d^4u}{ds^4} + \kappa^2\right]u = 0$$

then I can obtain an approximation using perturbation methods, i.e., the eigenvalues are approximated by

$$\kappa_n^2 = n^2\pi^2 + \epsilon n^4\pi^4$$

for small $\epsilon$. Then, for small values of $\epsilon$, I can use the values $\kappa_n^2$ as initial guesses to the Newton algorithm. But when $\epsilon$ increases these initial guesses fail.

Since I know the original differential equation, I can use FEM or FDM to find the eigenvalues, but I am interested in other methods. Below, you can see the comparison using FDM (1001 points), perturbation solution and Newton method (using the perturbation solution as initial guess). All the curves are for $\epsilon=10^{-3}$, but I couldn't make those guesses to work for greater values of $\epsilon$.

enter image description here

Question: Is there any other method to solve this problem? Maybe some kind of transformation that can be applied to the equation?

$\endgroup$
6
$\begingroup$

I suspect the main problem is the magnitude of the values. If you divide through by $\cosh$ to make all the numbers smaller, then ApproxFun doesn't seem to have a problem finding all the roots.

function lhs(ε, α)
    A, B = sqrt((α-1)/(2ε)), sqrt((α+1)/(2ε))
    return sqrt(α*α-1)*(cos(A) - 1/cosh(B)) - sin(A)*tanh(B)
end

u = Fun(x -> lhs(1e-3, x), [1.0, 10.0], 1024)
r = ApproxFun.roots(u)
r2 = [fzero(x -> lhs(1e-3, x), r[i]*(1+1e-6*rand())) for i in 1:length(r)]

This finds all the roots in $\alpha\in[1,10]$, for $\epsilon=10^{-3}$, $\kappa^2\lesssim 25000$, using an 1024-term Chebyshev series, taking a very small amount of time. This seems to be accurate to about $10^{-8}$ and the roots can be polished off with Newton's method or similar.

$\endgroup$
  • $\begingroup$ Thanks for the answer. Is there a Python equivalent for ApproxFun? Indeed dividing by $\cosh$ make the numbers smaller, what translates in better solutions. Just by doing that I managed to move up to $\epsilon\sim 10^{-2}$. Also, solving for $\alpha$ instead of $\kappa^2$ let me go up to $\epsilon\sim 1$, so probably is wiser to solve for $\alpha$ instead. But, even after those changes I can go upper than that, maybe is because of the solver... $\endgroup$ – nicoguaro Jun 30 '16 at 1:38
  • 2
    $\begingroup$ @nicoguaro: The insight to Kirill's answer is that roots of univariate polynomials are easy to compute. If we can approximate any arbitrary univariate function as a univarite polynomial, then approximate roots are also easy to compute. Kirill performs the polynomial approximation numerically, but looking at your expression, it's not difficult to see that a polynomial approximation is also available analytically, e.g. using Padé approximants. $\endgroup$ – Richard Zhang Jun 30 '16 at 3:35
  • 2
    $\begingroup$ @nicoguaro According to chebfun.org/about/projects.html, there is pychebfun which seems like the same thing in python. $\endgroup$ – Kirill Jun 30 '16 at 6:14
  • $\begingroup$ @Kirill, I tried with pychebfun and it works nicely. One more time, it seems a better idea to solve for $\alpha$ and not for $\kappa^2$. $\endgroup$ – nicoguaro Jun 30 '16 at 15:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.