enter image description hereenter image description hereSo I am trying to recreate the numerical treatment in this papir (https://journals.aps.org/prb/abstract/10.1103/PhysRevB.29.130), equation 3.7. The main point is to try to solve the following equation

$$z(\omega)=\frac{1}{1+2*\alpha|\omega|+\omega^2}*\lambda*\int_{-\infty}^{\infty} d\omega' z(\omega-\omega')z(\omega'),$$

where $\lambda$ is a variable we define through and $\alpha$ is a constant we choose. There is a trivial solution to this which is $z(\omega)=0$, which I want to avoid. We know that $z(\pm\infty)=0$ and that it is symmetric around the y-axis. The procedure in the article is an iteration procedure

  1. Start with an initial $z_0(\omega)$ and $\lambda_0=\frac{3}{2\sqrt{2}\Pi}$
  2. Calculate $z_1(\omega)=O(\lambda_0,z_0(\omega))$, where O is the operator defined by the right hand side of the equation
  3. calculate $\lambda_1 = \lambda_0\zeta^{-2}$, where $\zeta = \frac{z_1(\omega=0)}{z_0(\omega=0)}$
  4. repeat, thereafter find $z_2(\omega)=O(\lambda_1,z_1(\omega))$ and so on. untill we have a converged function.

I am trying to perform this for $\alpha=1$, which should give us an shape $\propto exp(-\omega)$. My initial code is below for one iteration:

        approxsoln [x_] = 1;
        \[Lambda] = 3/(2*Sqrt[2]*Pi);
         values = Table[{x, 1/(1 + 2*1*x + x^2)*\[Lambda]
          NIntegrate[approxsoln[y] approxsoln[x - y],
             {y, -Infinity, Infinity}]}, {x,0, 2, 1/100}];
        approxsoln1[x_] = InterpolatingPolynomial[values, x];
        \[Lambda]1 = \[Lambda]*(approxsoln1[0]/approxsoln[0])^(-2)

And then I will just repeat my code. My problem is to deal with the $|\omega|$ and automating it.

Hope someone can help


I have uploaded a picture of the solution they got, which is normalised


3 Answers 3


Here is one way (using NestWhile or NestWhileList):

Define my "operator" o. Of course you should modify this to make sure it meets your convergence needs (such as increasing the range and number of sampling points).

     zn = Interpolation[
         {#,1/(1+2 α Abs[#]+ #^2)λ NIntegrate[z[#-w]z[w],{w,0,2}]}&/@Range[0,2,1/20]
     {λ (z[0]/zn[0])^2,zn}

You could choose another test of course.

res = NestWhileList[o @@ # &, {3/(2 Sqrt[2] Pi), Function[x, 1]}, 
  Echo[NIntegrate[Abs[#[[2]][x] - #2[[2]][x]], {x, 0, 2}]/
      NIntegrate[#2[[2]][x], {x, 0, 2}]] >  10^-2 &, 2, 20]

  Evaluate@Through[res[[{1, 2, 6, 11, 16, 21}, 2]][w]], {w, 0, 0.5},
  PlotLegends -> StringTemplate["n=``"] /@ ToString /@ ({1, 2, 6, 11, 16, 21} - 1),
  PlotRange -> All]

Mathematica graphics


here is a loop approach. This can be fixed up to use NestWhile, but I think it easier to see what happens here..

Note we look at a finite domain due to use of interpolating functions (there is probably a way around that.. )

rmax = 1000
o[a_, l_, z_, x_?NumericQ] := 
 1/(1 + 2 a Abs[x] + x^2) l NIntegrate[
   z[y] z[x - y], {y, -rmax, rmax}]

lam = 3/(2 Sqrt[2] Pi);
f0 = Exp[-Abs[#]] &;  (* note this is a better guess but f0=1& works ok*)
Monitor[ff = Table[
    f1 = FunctionInterpolation[
           o[ 1  , lam , f0 , x ], {x, -rmax, rmax}];
    lam = lam/(f1[0]/f0[0])^2;
    f0 = f1, {i, 20}], {i, lam}];

Plot[ff[[-1]][x], {x, -1000, 1000}, PlotRange -> All]

enter image description here

  • $\begingroup$ I do not like that it goes to minus, I have uploaded a picture showing what they got for $\alpha=1$. here it seems that every solution for every alpha is the same $\endgroup$ Commented Apr 2, 2018 at 21:06
  • $\begingroup$ well, I think I captured the algorithm you spelled out, not sure whats wrong ( no permissions to see that paper.. ) $\endgroup$
    – george2079
    Commented Apr 2, 2018 at 21:34
  • $\begingroup$ I do not know if it helps, but I have attach the relevant part of the article, if you can see a mistake with my method $\endgroup$ Commented Apr 2, 2018 at 21:42
  • $\begingroup$ I checked your solution out of curiosity and found that FunctionInterpolation was the cause of the problem. For the time being, let me create another answer for better communication. Since it is a mere edit of your post, please feel free to change and incorporate into your answer. Eventually, I can delete my post. $\endgroup$
    – Sungmin
    Commented Apr 2, 2018 at 22:20
  • $\begingroup$ @sungmin no problem, good call. $\endgroup$
    – george2079
    Commented Apr 3, 2018 at 13:37

This is a temporary answer to give some feedback to george2079's answer.

It seems that the root cause of the problem in his post is in the use of FunctionInterpolation. Since I do not know the implementation detail of the Function, I cannot comment on why it happens.

To see the problem, check the following code snippets,

makeInter[ lam_, f0_] := 
 FunctionInterpolation[o[1, lam, f0, x], {x, -rmax, rmax}]
makeInter2[ lam_, f0_] := 
 Table[{x, o[1, lam, f0, x]}, {x, -rmax, rmax, 0.01}]

tmp = makeInter[lam, f0] // Quiet
data = makeInter2[lam, f0];

Plot[{xx[x], Interpolation[data][x]}, {x, -5, 5}, PlotRange -> All]

The only difference between two functions makeInter and makeInter2 is the explicit generation of data used for a interpolation. One can check that they do not produce the same curves especially when rmax is very large.

By the way, here is the working code (I set rmax =5 which turns out to be enough for the present goal.) It seems to be quite close to the original plot in the paper.

rmax = 5;
o[a_, l_, z_, x_?NumericQ] := 
 1/(1 + 2 a Abs[x] + x^2) l NIntegrate[z[y] z[x - y], {y, -rmax, 

lam = 3/(2 Sqrt[2] Pi);
f0 = Exp[-Abs[#]] &;(*note this is a better guess but f0=1& works ok*)

Monitor[ff = Table[f1 = Interpolation[makeInter2[lam, f0]];
    lam = lam/(f1[0]/f0[0])^2;
    f0 = f1, {i, 20}], {i, lam}];

lamU = 3/2 1/Sqrt[2 \[Pi]];
Plot[{lam/lamU ff[[-1]][x]}, {x, 0, 2}, PlotRange -> {0, 5}]

enter image description here

  • $\begingroup$ is it normal to get a lot of red warnings in the beginning? $\endgroup$ Commented Apr 2, 2018 at 22:40
  • $\begingroup$ @MortenSode It is an expected behavior because the interpolated range is finite but the integral is performed over the entire real line. Since the function vanishes quickly, this should not change the result significantly. By the way, the warnings can be ignored by adding Quiet in the function o. $\endgroup$
    – Sungmin
    Commented Apr 2, 2018 at 23:02
  • $\begingroup$ perfect thank you... I will try to run it, but it seems to take a long time...How long did it take for you? $\endgroup$ Commented Apr 2, 2018 at 23:06
  • $\begingroup$ @MortenSode It takes several minutes on my machine. Since chuy's answer controls the number of iterations, it could be faster than this. $\endgroup$
    – Sungmin
    Commented Apr 2, 2018 at 23:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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