iteratively solve integral equation

Updated:

I am solving the following integral equation for $f(k)$

$$f(k)=-\iint d^2p\,d^2q\frac1{g(p)}\frac1{g(p-q)}\left(\frac 1{g(k-q)}-\frac1{g(-q)}\right)$$

where $g(k)=k^2+f(k)$

so the equation is

$$\small f(k)=-\iint d^2p\,d^2q\frac1{p^2+f(p)}\frac1{(p-q)^2+f(\vert p-q\vert)}\left(\frac 1{(k-q)^2+f(\vert k-q\vert)}-\frac1{q^2+f(q)}\right)$$

where $f(k)$ and $g(k)$ are isotropic functions in 2D.

Before I start numerically solving it, I am expecting that $f(k)$ is linearly increasing at small $k$ from $(0,0)$ and becomes constant for large $k$.

Here I tried to solve it iteratively, starting with the trial function $f(k)=1$:

f[k_] = 1;
g[k_] = k^2 + f[k];
iterstep := (values = Table[{k, NIntegrate[-p q/
g[p]/(p^2 + q^2 - 2 p q Cos[ϕp - ϕq] + f[Sqrt[ p^2 + q^2 - 2 p
q Cos[ϕp - ϕq] ]]) (1/(k^2 + q^2 - 2 k q Cos[ϕq] +
f[Sqrt[k^2 + q^2 - 2 k q Cos[ϕq]]])-1/g[q]), {p, 0, 50}, {q, 0,
50}, {ϕp, 0, 2 π}, {ϕq, 0, 2 π}, Method -> "QuasiMonteCarlo",
PrecisionGoal -> 4,]}, {k, 0, 50, 10}] ;
f1[k_]= InterpolatingPolynomial[values, k];
f[x_] = Piecewise[{{f1[x], x < 50}, {f1, x > 50}}]
g[k_] = k^2 + f[k];)
plot := Show[Plot[f[k], {k, 0, 50}, PlotRange -> All], ListPlot[values]];

I do the integral from 0 to a cutoff 50 and evaluate $f(k)$ at 5 points form 0 to 50 then do a fit to get new function $f(k)$ for next step.

after 1st step:

iterstep
values
plot 2nd step 3rd step 4th step The major problem now: iteration is not contractive. Jump between high slope to low slope to even higher and even lower. How to change the program? Another better method?

The minor problem: I want to improve numerical integral accuracy. I have been asking about related numerical integral before: Multidimensional NIntegral with singularity. Some good advices were given.

This is more a couple of comments than an answer, but it requires some space. So, I post it here.

First: I observed that you use InterpolationPolynomial. Due to Runge's phenomenon, this is a really bad idea when interpolating so many data points. Better use Interpolation instead. This uses piecewise polynomials.

Second: Your function definitions are a bit odd (no Blanks, no SetDelayed). I took the freedom to change that (see below).

Third: What are you actually doing there? It looks to me like a fixed point iteration that simply reevaluates the right hand side of the equation again and again. Do you have evidence that your iteration operator is contractive? Being noncontractive might be the reason why this does not work...

f = 1. &;
g[k_] := k^2 + f[k];
iterstep[] := (values =
Table[
Print[k]; {k,
NIntegrate[-p q/
g[p]/(p^2 + q^2 - 2 p q Cos[\[Phi]p - \[Phi]q] +
f[Sqrt[p^2 + q^2 - 2 p q Cos[\[Phi]p - \[Phi]q]]]) (1/(k^2 +
q^2 - 2 k q Cos[\[Phi]q] +
f[Sqrt[k^2 + q^2 - 2 k q Cos[\[Phi]q]]]) - 1/g[q]),
{p, 0, 50},
{q, 0, 50},
{\[Phi]p, 0, 2 \[Pi]},
{\[Phi]q, 0, 2 \[Pi]}
]
},
{k, 0, 75, 5}];
f = Interpolation[values];
)
• thanks for the reply. what is "[ ]" for after iterstep? About your third comment, your understanding on my code is right. whether it is contractive or not, I am not sure yet. I think I have not make the Nintergral done properly.
– p.s
Jul 31 '17 at 23:17
• @p.s: The brackets [] thurn iterstep into a function without arguments. I did that just for a habit. This prevents the large computation to be started, when only the name of the function gets evaluated. Aug 1 '17 at 6:20
• @ Henrik Schumacher I updated post. Looks like the iteration is not contractive.
– p.s
Aug 2 '17 at 2:02