0
$\begingroup$

Possible Duplicate:
Kramers-Kronig in Mathematica

I am trying to get epsilon2 (imaginary part) from known epsilon1 (real part) data with the respective energy values (w), using the Kramers-Kronig relations. My code is:

{epsilon1, w} = ToExpression@Import["http://pastebin.com/raw.php?i=Z08GwMa8"];
epsilon2[w_]:= 2/Pi NIntegrate[a epsilon1[a]/(a^2-w^2), {a, 0, 5.5793}, 
    Method-> "PrincipalValue", MaxRecursion->20, Exclusions->{(a^2-w^2)==0}]//Quiet

Plot[epsilon2[w],{w,0,5.57932},AxesOrigin->{0,0},PlotPoints->400]

However, by using the code above, I managed only to get a list of irrelevant numbers. The graph of epsilon1 looks like a sine curve, while that of epsilon2 looks like an "M" shape. Any help is greatly appreciated! =)

P.S. If pastebin is down, you can get the original data from revision 1.

Improve this question
$\endgroup$
0

1 Answer 1

Reset to default
4
$\begingroup$

epsilon1as you defined it, isn't a function but a list.
You may try something like:

{epsilon1, w} = ToExpression@Import["http://pastebin.com/raw.php?i=Z08GwMa8"];
f = Interpolation[Transpose[{Flatten[w], Flatten[epsilon1]}]]
epsilon2[w_] := epsilon2[w] = 
                2/Pi NIntegrate[a f[a]/(a^2 - w^2), {a, 0, 5.5793}, 
                Method -> "PrincipalValue", Exclusions -> {(a^2 - w^2) == 0}] //  Quiet
DiscretePlot[epsilon2[w], {w, 0.5, 5.5, .25}, AxesOrigin -> {0, 0}, Joined -> True]

Mathematica graphics

Improve this answer
$\endgroup$
4
  • $\begingroup$ Hi @belisarius! Thank you =) 2 things: 1) The shape of epsilon2 is there, but it is supposed to be above the x-axis (i.e. flipping the negative region up to the positive region). 2) From my epsilon2 data, the curve should start off from the origin {0,0} rather than from the high values. Is there a problem or something? =) $\endgroup$
    – weesiang
    Oct 16, 2012 at 2:27
  • $\begingroup$ @weesiang You may have to check them. That curve is a reasonable approach to your integrand. $\endgroup$
    – Dr. belisarius
    Oct 16, 2012 at 10:22
  • 1
    $\begingroup$ Oh yes!! I miss out a parameter in the equation! Now I got it! Thank you very much @belisarius! $\endgroup$
    – weesiang
    Oct 16, 2012 at 10:47
  • $\begingroup$ Hi @belisarius, If I were to integrate from 0 to infinity, how would you alter the syntax? And also, I would like to ask you if there are alternative commands to DiscretePlot? =) $\endgroup$
    – weesiang
    Dec 6, 2012 at 11:20

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