# integral equation iteration

So 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

Edit:

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

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).

o[λ_,z_]:=
Block[{α=1,zn},
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]

Plot[
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]


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]


• 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 – Morten Sode Apr 2 '18 at 21:06
• well, I think I captured the algorithm you spelled out, not sure whats wrong ( no permissions to see that paper.. ) – george2079 Apr 2 '18 at 21:34
• 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 – Morten Sode Apr 2 '18 at 21:42
• 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. – Sungmin Apr 2 '18 at 22:20
• @sungmin no problem, good call. – george2079 Apr 3 '18 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,
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}]


• is it normal to get a lot of red warnings in the beginning? – Morten Sode Apr 2 '18 at 22:40
• @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. – Sungmin Apr 2 '18 at 23:02
• perfect thank you... I will try to run it, but it seems to take a long time...How long did it take for you? – Morten Sode Apr 2 '18 at 23:06
• @MortenSode It takes several minutes on my machine. Since chuy's answer controls the number of iterations, it could be faster than this. – Sungmin Apr 2 '18 at 23:09