13

I am trying to calculate the change of the refractive index from the change of the absorption coefficient using the Kramers-Kronig relations, in Mathematica.

c = 300000000;

daF[l_] = 500 * 0.28 Exp[-((l - 500)/90)^2];

dnFpoints = Table[
    {
        ln,
        c/Pi NIntegrate[
            daF[li] / ((2 Pi c 10^9 /li)^2 - (2 Pi c 10^9 / ln)^2),
            {li, 800, 200},
            Method -> {"PrincipalValue"},
            Exclusions -> ((2 Pi c 10^9 /li)^2 - (2 Pi c 10^9 / ln)^2) == 0
        ]
    },
    {ln, 300, 600}
];

Unfortunately, Mathematica displays an error that it does not converge to prescribed accuracy and the output is junk (I would expect a smooth curve with a negative minimum first and then a positive maximum). I am using version 8, if it matters. Any ideas?

Artes
  • 57,212
  • 12
  • 157
  • 245
mcandril
  • 633
  • 4
  • 20

1 Answers1

13

I think you intended to use {li, 200, 800} instead of {li, 800, 200}.

If you do so, then you could visualize the result :

ListLinePlot@dnFpoints

enter image description here

Moreover I would rather define daF in the following form :

daF[l_]:= 500 * 0.28 Exp[-((l - 500)/90)^2]
c = 3 10^8;

Edit

Instead of using Table of dnFpoints I add an alternative method for calculation of dnF function.

dnF[ln_] := 
  1/( 4c Pi^3 10^18 ) NIntegrate[ daF[li] / ( 1/li^2 - 1/ln^2 ), 
                                  { li, -\[Infinity],  ln, \[Infinity] }, 
                                  Method ->  "PrincipalValue", 
                                  Exclusions ->  Automatic  
                                ] // Quiet

In general one should choose appropriate options for NIntegrate like PrecisionGoal and MaxRecursion however in this case it is quite sufficient to use Quiet for evaluating of dnF function without outputting any messages generated.

Now we can plot dnF function increasing appropriately a range of the dependent variable, e.g. :

Plot[dnF[ln], {ln, 30, 900}, AxesOrigin -> {0, 0}, PlotPoints -> 200]

enter image description here

Artes
  • 57,212
  • 12
  • 157
  • 245
  • When you change ln for li what is ln then? I couldn't find a definition for it in the above code and get a non-numerical values error. – Matariki Feb 14 '12 at 18:43
  • 1
    @Matariki mcandril integrates numerically daF[li]/(f1(li)-f2(ln)) over li from 200 to 800, and then collects results in Table[{ln, integral(ln)}, {ln, 300, 600}]. Therefore we can't change variables li and ln. – Artes Feb 14 '12 at 18:57
  • Doh, missed that! Thanks for clarifying. – Matariki Feb 14 '12 at 19:08
  • ...if the integral he needs actually does go from 800 to 200, then it's a simple matter of changing signs, since $\int_a^b=-\int_b^a$. – J. M.'s missing motivation Feb 14 '12 at 23:21
  • 2
    Wow, thanks. I indeed meant 800-200, this comes from the substitution $\omega$->$\lambda$. Of course I knew J. M. comment, but I thought Mathematica does, too ... Is that a bug, or is there more behind the problem? – mcandril Feb 15 '12 at 07:58
  • @Artes Can you shortly explain why there has to be an error? I think the integral is defined in both ways, as J. M. explained. Also though I don't have performance issues, I do have issues with the result. I differs by 13 orders of magnitude from the MathLab- and the experimental value. I would expect the extrema in an order of magnitude between 10^-4 and 10^-6. As the schape of curve itself looks fine though, I don't think this is an issue of the numerics - though I cannot find an error in my formula. – mcandril Feb 16 '12 at 12:59
  • @mcandril If you evaluate e.g. NIntegrate[1/(x - x^2), {x, 2, -1}, Method -> "PrincipalValue", Exclusions -> x - x^2 == 0] the result is the same (modulo -) as in the case {x, -1, 2} unlike for the integrand in your case. I suppose the problem arises when numerically integrating in a "reverse direction" because of options Method -> "PrincipalValue", Exclusions -> something for certain kinds of functions. It is just a kind of issue when working with numerics, and do not think it is a bug since you get messages. – Artes Feb 16 '12 at 21:06
  • @mcandril I don't know why you put 10^9 in the integrand so I think this is the reason of appearing different orders of magnitude. – Artes Feb 16 '12 at 21:06
  • @Artes Thanks for the explanation. I work on the optical nm-scale. The Problem was that I forgot thet $\omega\to\lambda$ sustitution for $\text{d}\omega$. Awkward. – mcandril Feb 20 '12 at 10:30