# Solution of an implicit integral equation

The present question is pretty similar to this one.

Here I want to find the solution of the following implicit integral equation: $$P_H(u)=\int_{-\infty}^\infty dJ \int_{-\infty}^\infty du_1 \int_{-\infty}^\infty du_2 \frac{e^{-J^2}}{\sqrt{\pi}}P_H(u_1)P_H(u_2)\delta\bigg(u-\text{sign}\big(J(H+u_1+u_2)\big)\text{min}\big(|J|,|H+u_1+u_2|\big)\bigg)\,,$$ where $$P_H(u)$$ is parametrized by $$H\in\mathbb{R}$$ and it is normalized (it is a probability density function): $$\int_{-\infty}^\infty du P_H(u) = 1\,.$$ Following the linked well-answered question I tried to solve it iteratively, with the following code:

PB[H_?NumericQ] :=
With[{L = 20, g = 0.96},
NestList[
Function[fa,
Block[{int},
int[x_?NumericQ] :=
Block[{s},
NIntegrate[
E^-J^2/Sqrt[Pi]
DiracDelta[
x - Sign[J (H + u1 + u2)] Min[Abs[J],
Abs[H + u1 + u2]]] fa[u1]*fa[u2], {J, -L, L}, {u1, -L,
L}, {u2, -L, L},
Method -> {{Automatic, "LocalAdaptive"}[[-1]],
"SymbolicProcessing" -> 0}]];
Interpolation[
Table[{x, g fa[x] + (1 - g) int[x]}, {x, -L, L, .1}]]]],
E^-#^2/Sqrt[Pi] &, 5]];

(*Plot*)
H = 0.358;
approx = PB[H];
Show[{Plot[{Table[approx[[i]][x], {i, 1, 5}]}, {x, -10, 10}, PlotStyle -> {Black(*,{Dashed,Gray}*)}, PlotRange -> All]}]



The final output is a function which is not normalized, do you have any suggestion? Maybe use the cumulative function instead of the density itself, or maybe multiplying by the normalization at each step?

Edit I think one of the problems was the DiracDelta inside the NIntegrate. To solve this issue I used a UnitStep function to distinguish the cases $$|J|<|H+u_1+u_2|$$ and $$|J|>|H+u_1+u_2|$$ and then using the DiracDelta to solve one of the integrals. The resulting code is the following (where I also used ParallelTable):

PB[H_?NumericQ] :=
With[{L = 10, g = 0.5},
NestList[
Function[fa,
Block[{int},
int[x_?NumericQ] :=
Block[{s},
NIntegrate[
SetPrecision[
E^-x^2/Sqrt[Pi] fa[u1]*fa[u2]*
UnitStep[Abs[H + u1 + u2] - Abs[x]], 21], {u1, -L,
L}, {u2, -L, L}, WorkingPrecision -> 20] +
NIntegrate[
SetPrecision[
E^-J^2/Sqrt[Pi] fa[Sign[J] x - H - u2]*fa[u2]*UnitStep[Abs[J] - Abs[x]], 21], {u2, -L,
L}, {J, -L, L}, WorkingPrecision -> 20]];
Interpolation[
ParallelTable[{x, g fa[x] + (1 - g) int[x]}, {x, -L,
L, .1}]]]], E^-#^2/Sqrt[Pi] &, 4]];

(*Plot*)
H = 0.358;
approx = PB[H];
Show[{Plot[{ParallelTable[approx[[i]][x], {i, 1, 4}]}, {x, -3, 3},
PlotStyle -> {Black(*,{Dashed,Gray}*)}, PlotRange -> All]}]


I didn't implement the normalization, otherwise the needed time to get the output is too large. My problem now is indeed the code execution time, too large.

• It looks like variables J, u, u1, u2 and parameter H are reals. So maybe Exp[-u^2]/Sqrt[Pi] is exact solution? Commented Jun 22 at 8:48
• I don't think so, for two reasons: 1) setting as input a gaussian function, the output is always equal to the input whatever sigma and mean I use for the gaussian function (with the same problem of normalisation) 2) I tried to evaluate numerically the integral setting $P_H(u)=Exp(-u^2)/Sqrt(Pi)$ for a fixed $u$ but the result is not coherent Commented Jun 23 at 14:04
• Nevertheless, please, see my answer:) Commented Jun 24 at 16:20
• There are a lot symmetries that can be used for simplification. First of all, it is obvious that $P_H(u)=P_H(-u)$. Thus, you can get rid of sign inside the $\delta$ function. Commented Jun 24 at 22:16

## 2 Answers

To solve the implicit integral equation for ( P_H(u) ) iteratively, ensuring the solution is normalized at each step, you can follow a similar approach as in the linked question but include a normalization step after each iteration.

Here's an improved version of your code:

1. Normalization Function: normalize[fa_] computes the normalization factor and then normalizes the function ( fa ).
2. NestList for Iterative Approach: The NestList function iteratively applies the transformation to fa, which is initially set to the Gaussian ( e^{-x^2}/\sqrt{\pi} ).
3. Normalization at Each Step: The normalize function ensures that ( P_H(u) ) remains a probability density function at each iteration by dividing by the integral of ( fa ).
4. Plot: The plot shows the evolution of the approximation over iterations, with legends indicating each iteration.

This method should help you obtain a normalized solution for ( P_H(u) ) iteratively. You can adjust the number of iterations and other parameters as needed.

Here's an improved version of your code:

ClearAll[PB]

PB[H_?NumericQ] :=
Module[{L = 20, g = 0.96, int, normalize},
normalize[fa_] :=
With[{norm = NIntegrate[fa[x], {x, -Infinity, Infinity}]},
Function[x, fa[x]/norm]
];
NestList[
Function[fa,
Block[{int},
int[x_?NumericQ] :=
Block[{s},
NIntegrate[
Exp[-J^2]/Sqrt[Pi]
DiracDelta[
x - Sign[J (H + u1 + u2)] Min[Abs[J],
Abs[H + u1 + u2]]] fa[u1]*fa[u2], {J, -L, L}, {u1, -L,
L}, {u2, -L, L},
Method -> {"LocalAdaptive",
"SymbolicProcessing" -> 0}]];
normalize[
Interpolation[
Table[{x, g fa[x] + (1 - g) int[x]}, {x, -L, L, .1}]
]]
]
],
Exp[-#^2]/Sqrt[Pi] &, 5
]
]

(* Plot *)
H = 0.358;
approx = PB[H];
Show[
Plot[{approx[[#]][x] & /@ Range[Length[approx]]}, {x, -10, 10},
PlotStyle -> Table[ColorData[1, "ColorList"][[i]], {i, 1, 5}],
PlotRange -> All,
PlotLegends ->
LineLegend[
Table["Approximation (r=" <> ToString[i] <> ")", {i, 1, 5}]
]
],
PlotRange -> All, PlotLabel -> "Approximations"]


• Thanks for your answer! However it seems like it doesn't work, since the output is basically the same as the input, whatever gaussian function you put as input. But it can be a step towards the solution, thank you. Commented Jun 21 at 15:46

Note, that on every step we have

NIntegrate[g fa[x]+(1-g) int[x],{x,-L,L}]=1


From this equation we can compute g as follows

g=(1- NIntegrate[int[x],{x,-L,L}])/(NIntegrate[fa[x],{x,-L,L}]-NIntegrate[int[x],{x,-L,L}])


Using this equation, we can modifier last code as follows

Get["NumericalDifferentialEquationAnalysis"];

pw = GaussianQuadratureWeights[50, -5, 5];

PB[H_?NumericQ] :=
With[{L = 5},
NestList[
Function[fa,
Block[{int, norm1, norm2, ff, ii},
int[x_?NumericQ] :=
Block[{u1, u2, J},
NIntegrate[
E^-x^2/Sqrt[Pi]  fa[u1]*fa[u2]*
UnitStep[(H + u1 + u2)^2 - x^2], {u1, -L, L}, {u2, -L, L},
Method -> "AdaptiveQuasiMonteCarlo", AccuracyGoal -> 2,
PrecisionGoal -> 2] +
NIntegrate[
E^-J^2/Sqrt[Pi]  fa[Sign[J]  x - H - u2]*fa[u2]*
UnitStep[J^2 - x^2], {u2, -L, L}, {J, -L, L},
Method -> "AdaptiveQuasiMonteCarlo", AccuracyGoal -> 2,
PrecisionGoal -> 2]]; ff = fa[#[[1]]] & /@ pw;
ii = int[#[[1]]] & /@ pw; norm1 = pw[[All, 2]] . ff;
norm2 = pw[[All, 2]] . ii; g = (1 - norm2)/(norm1 - norm2);
Interpolation[Transpose[{pw[[All, 1]], g  ff + (1 - g)  ii}]]]],
E^-#^2/Sqrt[Pi] &, 4]];


We used Gauss quadrature instead of NIntegrate. Visualization

H = 0.358;
approx = PB[H]; // AbsoluteTiming

(*Out[]= {5.85601, Null}*)

Show[Plot[E^-x^2/Sqrt[Pi], {x, -3, 3}, PlotStyle -> Blue],
Plot[Evaluate[Table[approx[[i]][x], {i, 2, 5}]], {x, -3, 3},
PlotRange -> {0, .6}, Frame -> True, PlotLegends -> Automatic,
PlotStyle -> Dashed]]
`

We see that initial condition is probably exact solution. And we can try to prove it. First, we note that same result computed by zeraoulia rafik as well.