Sign up ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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[""];
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


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.

share|improve this question

marked as duplicate by whuber, Oleksandr R., rcollyer, R. M. Oct 27 '12 at 3:17

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

1 Answer 1

up vote 3 down vote accepted

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

{epsilon1, w} = ToExpression@Import[""];
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

share|improve this answer
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? =) – weesiang Oct 16 '12 at 2:27
@weesiang You may have to check them. That curve is a reasonable approach to your integrand. – belisarius has settled Oct 16 '12 at 10:22
Oh yes!! I miss out a parameter in the equation! Now I got it! Thank you very much @belisarius! – weesiang Oct 16 '12 at 10:47
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? =) – weesiang Dec 6 '12 at 11:20

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