0
$\begingroup$

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$):

https://i.sstatic.net/ZV8Rs.png

Improve this question
$\endgroup$
7
  • $\begingroup$ Define a function in terms of alpha, as here. $\endgroup$ Commented Jun 4, 2014 at 18:36
  • $\begingroup$ 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. $\endgroup$ Commented Jun 4, 2014 at 18:41
  • 1
    $\begingroup$ Personally, I would be more inclined to try to help if there were Mathematica code for the matrix I could copy and work with. $\endgroup$ Commented Jun 5, 2014 at 1:41
  • $\begingroup$ @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. $\endgroup$ Commented Jun 5, 2014 at 6:13
  • $\begingroup$ @george2079 $\beta$ is a number, $u$ is a function of the main variable $w$. I posted a picture of the code. $\endgroup$ Commented Jun 5, 2014 at 6:16

1 Answer 1

Reset to default
3
$\begingroup$

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 == D[Det[mat] /. {u -> u[α]}, α], (* DE -- see Note below *) 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}]; Show[ ParametricPlot[ Through[{Re, Im}[sol[α]]], {α, 1, maxalpha}, ColorFunction -> "Rainbow", AspectRatio -> 1/2, PlotRange -> All], Graphics[ {Point@roots} ] ]

Mathematica graphics

The precision of NDSolve (compared with FindRoot):

ListLogPlot[ 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 /. FindRoot[Det[mat] /. α -> 1, {u, 1}, WorkingPrecision -> wp]), α'[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:

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

Mathematica graphics

Note: Stelios pointed out that the original DE 0 == Dt@Det[mat] /. {u -> u[α]}, which works in V9/10, does not work in V8. So I updated it, so that the answer might work in other, hopefully all, versions of Mathematica.

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ Very elegant approach. Shouldn't the (*DE*) equation of the first method be 0 == D[Det[mat] /. {u -> u[\[Alpha]]}, \[Alpha]] though? $\endgroup$ Commented Jun 27, 2015 at 9:48
  • $\begingroup$ @Stelios I can't find any proof/documentation, but I believe NDSolve treats Dt[α] as 1, since α is the independent variable. So they're equivalent. An alternative explanation is that the Dt[α] cancels out when NDSolve solves for the derivative. In any case, your suggestion is perhaps the more normal way. (In PDEs, on the other hand, there's a big difference between Dt and partial derivatives. Dt might never be appropriate or work.) $\endgroup$ Commented Jun 27, 2015 at 14:18
  • 2
    $\begingroup$ Thank you for your detailed response. The reason I posted the comment was that my MMA (version 8.0) gives an error message with the Dt version of the equation, whereas it runs fine with the D[#,a] version, and wanted to make sure that the latter is indeed equivalent with your approach. $\endgroup$ Commented Jun 27, 2015 at 15:48
  • 1
    $\begingroup$ @Stelios Yes, I get the same thing both ways on V10.1 Given that it doesn't work in earlier version, I might ought to change this answer. (But I've upvoted your comments to help draw attention to the issue.) $\endgroup$ Commented Jun 27, 2015 at 16:05

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.