Take the 2-minute tour ×

I'm trying to derive some of the results of the following paper:

Electrodynamics of semiconductor-coated noble metal nanoshells, JT Manassah - Physical Review A

In the paper there is matrix $\mathbf M_l^M$, defined as follows, with $\alpha$ as a parameter:

enter image description here

as you know, $j$ and $n$ and $h$ are different types of spherical Bessel functions.

Now, I want t solve the following equation: $$\det ({\bf M}_l^M)=0$$

and obtain the roots in terms of $\alpha$.( I think the equation doesn't have analytical solution, and so I tried to use FindRoot command, but this equation has $\alpha$ as a parameter and the roots should be founded with respect to $\alpha$. If there is no (better) way other than FindRoot, how can I do this using this command? I mean, how can I make FindRoot to solve the equation for different $\alpha$s and what is $x_0$ in the FindRoot? (The roots are complex numbers.)

(All I want is to obtain a plot, showing the real part of the root versus $\alpha$)

Here is my code (definition of $\bf M$):


share|improve this question
Define a function in terms of alpha, as here. –  Szabolcs Jun 4 '14 at 18:36
you will probably do well to use the solution at each alpha as the starting point for the next point. Your expression appears to have a beta parameter as well.. by the way. –  george2079 Jun 4 '14 at 18:41
Personally, I would be more inclined to try to help if there were Mathematica code for the matrix I could copy and work with. –  Michael E2 Jun 5 '14 at 1:41
@MichaelE2 Yeah, I know; but the letters used in the code have subscripts, superscripts, radicals, etc., so I can't copy-paste the code. I posted a picture of it for now. If there is a better way for posting the code let me know. –  user215721 Jun 5 '14 at 6:13
@george2079 $\beta$ is a number, $u$ is a function of the main variable $w$. I posted a picture of the code. –  user215721 Jun 5 '14 at 6:16

1 Answer 1

I don't want to type in the OP's code and have to worry about typos, etc. But here's an example that in theory should show how to work with the OP's example. One can use NDSolve in one of two ways to get one variable u in terms of another α.

One issue that remains is the the function are oscillatory and the equation probably has many roots for a given α.

Here's my matrix, more or less randomly made of Bessel functions, etc:

mat = {
   {SphericalBesselJ[1, α u], -SphericalHankelH1[3, α u],  0},
   {SphericalBesselY[1, α u],  SphericalBesselJ[1, α u],  -SphericalHankelH2[3, α u]},
   {0,                         SphericalBesselY[1, α u],   SphericalBesselJ[3, α u]}

First Method

We can turn the equation into a differential equation by differentiating it and specifying an initial value (found with FindRoot). The precision wp may be set to MachinePrecision instead of 20. It will be faster, but the second method does not work with machine precision. The setting wp = 20 is used for the sake of comparison.

wp = 20;
maxalpha = 20;
sol = u /. First@NDSolve[{
      0 == Dt@Det[mat] /. {u -> u[α]},            (* DE *)
      u[1] == (u /.                               (* IV *)
         FindRoot[Det[mat] /. α -> 1, {u, 1}, WorkingPrecision -> wp])},
     {u}, {α, 1, maxalpha}, WorkingPrecision -> wp];

Using FindRoot with high precision for the sake of comparison:

roots = Table[
   Through[{Re, Im}[u /. FindRoot[Det[mat], {u, sol[α]}, WorkingPrecision -> 100]]],
   {α, maxalpha}];

  Through[{Re, Im}[sol[α]]], {α, 1, maxalpha},
  ColorFunction -> "Rainbow", AspectRatio -> 1/2, PlotRange -> All],

Mathematica graphics

The precision of NDSolve (compared with FindRoot):

 Abs /@ ((Table[sol[α], {α, maxalpha}] - roots.{1, I}) / roots.{1, I})

Mathematica graphics

Second method

The second way is to set up a DAE to trace the solution. This is slower and less accurate on this example. I suspect the reason is that NDSolve is more interested in the accuracy of the solution α to the DE than the accuracy of the solution u determined by the constraint. So there's not much to recommend it here. However, it can be a nice way to construct an interpolating function of a quantity that is a function of parametrized point on an integral curve.

wp = 20;
maxalpha = 20;
solDAE = u /. First@NDSolve[{
     0 == Det[mat] /. {u -> u[t], α -> α[t]},
     u[1] == (u /. 
     α'[t] == 1, α[1] == 1},
    {u, α}, {t, 1, maxalpha},
    Method -> "StateSpace", WorkingPrecision -> wp]

Comparison of precision, which is about two orders of magnitude worse than the first method:

 Abs /@ ((Table[solDAE[α], {α, maxalpha}] - roots.{1, I}) / roots.{1, I})

Mathematica graphics

share|improve this answer

Your Answer


