# Solution of an integral equation

I want to numerically solve the following integral equation for $$A(x)$$ (where $$A: \mathbb{R}\mapsto\mathbb{R}^+$$): $$A(x)=\exp{\left(-\int_{-x}^{\infty} (x+y)^r A(y) dy\right)},$$

where $$r\in\mathbb{N}$$ (I am interested in $$r=1,2$$). Moreover I know that $$A(-\infty)=1$$ and $$A(+\infty)=0$$. What I expect is a unique solution also for the cases $$r=1$$ and $$r=2$$.

I tried something using FindRoot:

precision = 4;
Maxit = 300;
accuracy = Round[precision/2];
r = 1;
fun[x_?NumericQ, A_?NumericQ] :=
A - Exp[-NIntegrate[(x+y)^r A, {y, -x, Infinity},
WorkingPrecision -> precision, AccuracyGoal -> accuracy,
PrecisionGoal -> precision,
Method -> {Automatic, "SymbolicProcessing" -> 0}]];
Afun[x_?NumericQ] :=
FindRoot[fun[x, z], {z, 0.2}, MaxIterations -> Maxit][[1]][[2]]
data = Table[{x, Afun[x]}, {x, -10, 10, 0.5}]
A = Interpolation[DeleteDuplicates[data]];
Plot[A[x], {x, -10, 10}] Afun[x_?NumericQ] :=
FindRoot[fun[x, z], {z, 0.2}, MaxIterations -> Maxit][[1]][[2]]


but this cannot clearly give the solution, since the function $$A(x)$$ depends on the integral variable. Can anyone help me?

Moreover I am referring to this question https://math.stackexchange.com/questions/4728120/solution-of-a-simple-integral-equation/4728209#4728209 . Here the equation for $$r=0$$ is solved analytically, but for $$r=1$$ and $$r=2$$ it seems impossible to tackle analytically the problem, so that I tried numerically. Maybe it is easier to solve the corresponding differential equation (which can be found deriving both sides).

• Where does condition A[Infinity]==0 come from? Jul 6, 2023 at 14:29
• In the problem I am studying $1-A(x)$ is the cumulative distribution function of a density which has support on the real line, that's why the two conditions. Jul 6, 2023 at 14:31
• DSolve can solve some integral equations exactly in Mathematica, but unsurprisingly (given the non-linearity) it returns unevaluated on this one. Unfortunately, it does not appear that NDSolve can do anything with integral equations. Jul 6, 2023 at 15:13
• Perhaps, with A[x,r]:= A[x]... , knowing A[x,0]=1/(1+Exp[x]), the recursion D[Log[A[x,r]],r]==r Log[A[x,r-1] , r,1,2,... might give a solution??? Jul 7, 2023 at 7:40
• Another idea might be an iterative approach. Therefore it would be helpful if the definition range of x,y and the integration range might be transformed in a finite intervall... Jul 7, 2023 at 9:43

A small modification of @AlexTrounev 's great answer gives a unique fixedpoint without jumps for r=0,1,2!

It's only necessary to modify the iterationpart $$\text{fa}_{n+1}(x)=e^{-\text{int}\left(\text{fa}_n(x)\right)}$$ .

Here I assume a memory based iteration $$\text{fa}_{n+1}(x)=(1-g) e^{-\text{int}\left(\text{fa}_n(x)\right)}+g \text{fa}_n(x)$$ (memory parameter 0<g<1) which works quite well.

Playing around with memory parameter g it seems g=1/4 is a good choice!

nlg[r_?NumericQ] :=
With[{L = 20, g = 1/4},
NestList[
Function[fa,
Block[{int},
int[x_?NumericQ] :=
Block[{s},
NIntegrate[(x + s)^r fa[s], {s, -x, L},
"SymbolicProcessing" -> 0}]];
Interpolation[
Table[{x, g fa[x] + (1 - g) Exp[-int[x]]}, {x, -L,
L, .1}]]]], (1 - Tanh[#])/2 &, 20]];


case r==0

rr = 0;
scale = 100000;
approx = nlg[rr];
Show[{ Plot[  {approx[[-1]][x](*,1/(1+Exp[x])*)}   , {x, -20, 20} ,
PlotStyle -> {Black(*,{Dashed,Gray}*)}, PlotRange -> All,
PlotLegends ->
LineLegend[{"approximation(r=" <> ToString[rr] <> ")"}]],
Plot[scale ( approx[[-1]][x] - approx [[-2]][x] )  , {x, -20, 20} ,
PlotStyle -> Red, PlotRange -> All,
PlotLegends ->
LineLegend[{Style["error * " <> ToString[scale], Red] }]]},
PlotRange -> All, PlotLabel -> "r=" <> ToString[rr]]


This solution agrees very well with the known exact solution 1/(1+Exp[x])!

case r==1

rr = 1;
scale = 100000;
approx = nlg[rr];
Show[{ Plot[  {approx[[-1]][x](*,1/(1+Exp[x])*)}   , {x, -20, 20} ,
PlotStyle -> {Black(*,{Dashed,Gray}*)}, PlotRange -> All,
PlotLegends ->
LineLegend[{"approximation(r=" <> ToString[rr] <> ")"}]],
Plot[scale ( approx[[-1]][x] - approx [[-2]][x] )  , {x, -20, 20} ,
PlotStyle -> Red, PlotRange -> All,
PlotLegends ->
LineLegend[{Style["error * " <> ToString[scale], Red] }]]},
PlotRange -> All, PlotLabel -> "r=" <> ToString[rr]]


case r==2

rr = 2;
scale = 100000;
approx = nlg[rr];
Show[{ Plot[  {approx[[-1]][x](*,1/(1+Exp[x])*)}   , {x, -20, 20} ,
PlotStyle -> {Black(*,{Dashed,Gray}*)}, PlotRange -> All,
PlotLegends ->
LineLegend[{"approximation(r=" <> ToString[rr] <> ")"}]],
Plot[scale ( approx[[-1]][x] - approx [[-2]][x] )  , {x, -20, 20} ,
PlotStyle -> Red, PlotRange -> All,
PlotLegends ->
LineLegend[{Style["error * " <> ToString[scale], Red] }]]},
PlotRange -> All, PlotLabel -> "r=" <> ToString[rr]]


The error( shown in red) between successive iterations is very small (no jumps!)

Hope it helps!

• Nice! Plotting the "error" is possible to tune the parameter $g$, for example for the case $r=2$ it seems that $g=.3$ gives the best result. I also tried to change the input function to see whether the iterative solution is "stable" and it seems ok. Thanks for the interest! Jul 10, 2023 at 10:26
• Tuning of g is a nice numerical task. More essential is the fact that my approach gives a unique stable fixedpointsolution for different values of g! I think that 's the detailed answer you asked for Jul 11, 2023 at 8:00
• @UlrichNeumann This is nice answer (+1). It looks like the average of 2 solutions. The question is, should there be one solution or 2? For example, a quadratic equation also has 2 solutions. Jul 11, 2023 at 13:13
• @AlexTrounev Thanks! Good question, from an iterative solution I would expect only one solution. That's why I tried to remember the current approximation (weight g ) in my approach Jul 11, 2023 at 13:24
• @AlexTrounev I tried a linear combination of your two solutions too, which gives c1 A1[x]+c2 A2[x]== A1[x]^c1*A2[x]^c2 . No idea how to find optimal parameters c1,c2 yet. Jul 11, 2023 at 13:30

To solve this problem we can use iterative method described here. Unfortunately two conditions at $$x=\pm \infty$$ are not enough to define unique solution. Numerically we can define 2 solutions as follows

With[{L = 20, r = 1},
nl = NestList[
Function[{fa},
Block[{int},
int[x_?NumericQ] :=
Block[{s},
r NIntegrate[ (x + s)^(r - 1) fa[s], {s, -x, L},
"SymbolicProcessing" -> 0}]];
NDSolveValue[{A'[x] == A[x] (-int[x]), A[-L] == 1},
A, {x, -L, L}]]], 1/(1 + Exp[#]) &, 10]];


Visualization

With[{L = 20}, {Plot[
Evaluate[Table[nl[[i]][x], {i, 6, 10, 2}]], {x, -L/2, L/2},
PlotLegends -> Table[i, {i, 6, 10, 2}], Frame -> True,
Axes -> False, GridLines -> Automatic],
Plot[Evaluate[Table[nl[[i]][x], {i, 5, 9, 2}]], {x, -L/2, L/2},
PlotLegends -> Table[i, {i, 5, 9, 2}], Frame -> True,
Axes -> False, GridLines -> Automatic]}]


As we can see solution jumps between two states. Same picture we have for $$r=2$$

With[{L = 20, r = 2},
nl1 = NestList[
Function[{fa},
Block[{int},
int[x_?NumericQ] :=
Block[{s},
r NIntegrate[ (x + s)^(r - 1) fa[s], {s, -x, L},
"SymbolicProcessing" -> 0}]];
NDSolveValue[{A'[x] == A[x] (-int[x]), A[-L] == 1},
A, {x, -L, L}]]], (1 - Tanh[#])/2 &, 10]];

With[{L = 20}, {Plot[Evaluate[Table[nl1[[i]][x], {i, 6, 10, 2}]], {x, -L/2, L/2},
PlotLegends -> Table[i, {i, 6, 10, 2}], Frame -> True,
Axes -> False, GridLines -> Automatic],
Plot[Evaluate[Table[nl1[[i]][x], {i, 5, 9, 2}]], {x, -L/2, L/2},
PlotLegends -> Table[i, {i, 5, 9, 2}], Frame -> True, Axes -> False,
GridLines -> Automatic]}]


Note that for $$r=0$$ we have functional equation with two solutions as well

With[{L = 20, r = 0},
nl0 = NestList[Function[{fa}, Block[{int}, int[x_?NumericQ] := fa[-x];
NDSolveValue[{A'[x] == A[x] (-int[x]), A[-L] == 1},
A, {x, -L, L}]]], 1/(1 + Exp[2 #]) &, 10]];

With[{L = 20}, {Plot[
Evaluate[Table[nl0[[i]][x], {i, 6, 10, 2}]], {x, -L/2, L/2},
PlotLegends -> Table[i, {i, 6, 10, 2}], Frame -> True,
Axes -> False, GridLines -> Automatic],
Plot[Evaluate[Table[nl0[[i]][x], {i, 5, 9, 2}]], {x, -L/2, L/2},
PlotLegends -> Table[i, {i, 5, 9, 2}], Frame -> True,
Axes -> False, GridLines -> Automatic]}]


Mean last 2 solutions (solid blue line) we can compare with exact solution $$A=1/(1+e^x)$$ (red dashed line) as

With[{L = 20},
Show[Plot[.5 (nl0[[9]][x] + nl0[[10]][x]), {x, -L/2, L/2},
Frame -> True, Axes -> False, GridLines -> Automatic,
PlotStyle -> Blue],
Plot[1/(1 + Exp[x]), {x, -L, L}, PlotStyle -> {Red, Dashed}]]]


Update 1 We can solve equation $$A(x)=e^{-\int_{-x}^{\infty}(x+y)^rA(y)dy}$$ directly using iterative method as follows

With[{L = 20},
nl0 = NestWhileList[
Function[{fa},
Block[{int},
int[x_?NumericQ] :=
NIntegrate[fa[y], {y, -x, L},
"SymbolicProcessing" -> 0}];
Interpolation[
Table[{x, Exp[-int[x]]}, {x, -L, L, .1}]]]], (1 - Tanh[#])/2 &,
Abs[#1[0] - #2[0]] > 10^(-9) &, 2, 20]];

nmax = Length[nl0];
With[{L = 20},
Plot[Evaluate[Table[nl0[[i]][x], {i, nmax-1,nmax}]], {x, -L/2,
L/2}, PlotLegends -> Automatic, Frame -> True,
GridLines -> Automatic]]


Here we used NestWhileList with convergence criteria, but have two solutions as well.

Code for arbitrary r

With[{L = 20, r = 1},
nl1 = NestList[
Function[fa,
Block[{int},
int[x_?NumericQ] :=
Block[{s},
NIntegrate[ (x + s)^r fa[s], {s, -x, L},
"SymbolicProcessing" -> 0}]];
Interpolation[
Table[{x, Exp[-int[x]]}, {x, -L, L, .1}]]]], (1 - Tanh[#])/2 &,
20]];


Update 2. As some explanation for Ulrich result $$g=1/4$$ we can use convergence criteria automated as follows

Table[With[{L = 20},
nl0 = NestWhileList[
Function[{fa},
Block[{int},
int[x_?NumericQ] :=
NIntegrate[fa[y], {y, -x, L},
"SymbolicProcessing" -> 0}];
Interpolation[
Table[{x, Exp[-int[x]] (1 - g) + g fa[x]}, {x, -L, L, .1}],
InterpolationOrder -> 4]]], (1 - Tanh[#])/2 &,
Abs[#1[0] - #2[0]] > 10^(-3) &, 2, 20]]; {g,
Length[nl0]}, {g, .21, .29, .01}]
(*{{0.21, 12}, {0.22, 10}, {0.23, 10}, {0.24, 8}, {0.25, 8}, {0.26,
8}, {0.27, 6}, {0.28, 10}, {0.29, 11}}*)


We can also plot this list to see the optimal g

ListLinePlot[%]


As we can see the optimal g=0.27 means minimal number of iterations to pass criterium $$|f_n(0)-f_{n-1}(0)|<10^{-3}$$. The corresponding solution for r=0,1,2

Do[With[{L = 20, g = 0.27},
nl[r] = NestWhileList[
Function[{fa},
Block[{int},
int[x_?NumericQ] :=
NIntegrate[(x + y)^r fa[y], {y, -x, L},
"SymbolicProcessing" -> 0}];
Interpolation[
Table[{x, Quiet[Exp[-int[x]]] (1 - g) + g fa[x]}, {x, -L,
L, .1}], InterpolationOrder -> 4]]], (1 - Tanh[#])/2 &,
Abs[#1[0] - #2[0]] > 10^(-3) &, 2, 20]];, {r, 0, 2, 1}]
With[{L = 20},
Plot[Evaluate[Table[nl[r][[-1]][x], {r, 0, 2, 1}]], {x, -L/2, L/2},
PlotLegends -> {0, 1, 2}, Frame -> True, GridLines -> Automatic]]



• Yours second code dosen't work for r=2 ? I have no Plot ?I tried with MMA 13.3. Jul 8, 2023 at 8:39
• @MariuszIwaniuk Thank you. Don't forget to add parameter L=20; before plot. Anyway post has been updated using With[]. Jul 8, 2023 at 9:52
• Ok works fine now. Thanks Jul 8, 2023 at 9:54
• @AlexTrounev Nice answer!!! Does iterative solution also jump for case r==0 ? Jul 8, 2023 at 10:36
• @UlrichNeumann Thank you for your code and question. I add case for $r=0$. As you can see the iterative solution jumps in this case as well. Jul 8, 2023 at 12:09

## Initial solution idea

Your integral equation can be reduced to an interesting system of differential equations, at least in $$r=1$$ case. Since $$A(x)>0$$ we can write $$A(x)=e^{a(x)}$$. Therefore $$a(x)=-\int_{-x}^{\infty}(x+y)e^{a(y)}dy.$$ Differentiate:

 D[-Integrate[(x + y) E^a[y], {y, -x, Infinity}], {x, 2}]


$$a''(x)=-e^{a(-x)}.$$ Let us now denote $$b(x)=a(-x)$$. In principle, we should demand continuity at $$x=0$$ $$a(0)=b(0)=c,\\ a'(0)=-b'(0)=d.$$ My thanks to @UlrichNeumann for correcting the sign in my second equation.

However, I do not know how to realize this numerically. Therefore, I will demonstrate a very rough approximation. Notice, that the motivation is to show a different approach. I am fully aware of the fact that there are issues that are hard to solve in general.

The first simplification is to replace the semi-infinite interval $$(-\infty,0]$$ with $$[-l, 0]$$.

The second simplification is to specify only the boundary values of functions and their derivatives at $$x=-l$$. We know that, in $$a(-\infty)=0, a'(-\infty)=0$$ and $$b(-\infty)=-\infty$$ and $$b'(-\infty)=0$$. In the finite interval, we replace infinite values with finite ones, selected by tries and errors (poor man shooting method) in such a way that $$a(0)\approx b(0)$$.

l = 5;
{u, v} = NDSolveValue[
D[a[x], {x, 2}] == -Exp[b[x]] && D[b[x], {x, 2}] == -Exp[a[x]] &&
a[-l] == 0 && a'[-l] == 0 && b[-l] == -20 && b'[-l] == 6.46, {a, b},
{x, -l, 0}]

w[x_] := v[-x]
g1 = Plot[{Exp[u[x]]}, {x, -l, 0}, PlotRange -> All,
PlotTheme -> {"Web", "Scientific"}];
g2 = Plot[{, Exp[w[x]]}, {x, 0, l}, PlotRange -> All,
PlotTheme -> {"Web", "Scientific"}];
Show[{g1, g2}]


## Improved solution for $$r=1$$

The question of @UlrichNeumann inspired me to generate a fully automatic and very accurate solution.

l = 5;
{u, v} =
NDSolveValue[{D[a[x], {x, 2}] == -Exp[b[x]],
D[b[x], {x, 2}] == -Exp[a[x]], a[-l] == 0, a'[-l] == 0,
a[0] == b[0], -a'[0] == b'[0]}, {a, b}, {x, -l, 0}]

A[x_] := Piecewise[{{Exp[u[x]], x <= 0}, {Exp[v[-x]], x > 0}}]


## Improved solution for $$r=2$$

Here I use the same method for $$a(x)=-\int_{-x}^{\infty}(x+y)^2e^{a(y)}dy.$$ A differential equation can be obtained by differentiating this identity 3 times.

D[-Integrate[(x + y)^2 E^a[y], {y, -x, Infinity}], {x, 3}]


yielding $$a'''(x)=-2 e^{a(-x)}$$. Introducing as before $$b(x)=a(-x)$$ I obtain \begin{aligned} a'''(x)&=-2 e^{b(x)},\\ b'''(x)&=+2 e^{a(x)},\\ a^{(n)}(-l)&=0,\\ b^{(n)}(0)&=(-1)^n a^{(n)}(0), \end{aligned} for $$n=0,1,2$$.

l = 5;
{u, v} =
NDSolveValue[{D[a[x], {x, 3}] == -2 Exp[b[x]],
D[b[x], {x, 3}] == 2 Exp[a[x]], a[-l] == 0, a'[-l] == 0,
a''[-l] == 0, a[0] == b[0], -a'[0] == b'[0], a''[0] == b''[0]}, {a,
b}, {x, -l, 0}]


• Your first Differentiate covers only the case r==0. This case has a general solution 1/(1+Exp[x])! Jul 10, 2023 at 13:11
• Note that for the $r = 2$ case the equation simply (?) becomes $a'''(x) = -2 e^{a(-x)}$. Jul 11, 2023 at 14:45
• @yarchik Thanks for re-presenting your approach, very interesting. Two remarks: #1 I think your second continuity condition should be -a'[0]==b'[0] . #2 I would expect NDSolveValue[{D[a[x], {x, 2}] == -Exp[b[x]], D[b[x], {x, 2}] == -Exp[a[x]], a[-l] == 1 , b[-l] == 0 , a[0] == b[0], -a'[0] == b'[0]}, {a, b}, {x, -l, 0} ] for the correct system? I tried the modifications but result isn't ok Jul 11, 2023 at 16:41
• @Yarchik Thanks for your reply. Your improved solution is great. Now it's possible to confirm the iterative approach r>0 (continous) for some values of r Jul 12, 2023 at 8:04
• @yarchik This is nice solution (+1). Jul 14, 2023 at 7:16