8
$\begingroup$

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

$\endgroup$
19
  • 1
    $\begingroup$ Where does condition A[Infinity]==0 come from? $\endgroup$ Commented Jul 6, 2023 at 14:29
  • $\begingroup$ 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. $\endgroup$ Commented Jul 6, 2023 at 14:31
  • $\begingroup$ 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. $\endgroup$ Commented Jul 6, 2023 at 15:13
  • 1
    $\begingroup$ 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??? $\endgroup$ Commented Jul 7, 2023 at 7:40
  • 1
    $\begingroup$ 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... $\endgroup$ Commented Jul 7, 2023 at 9:43

3 Answers 3

4
+50
$\begingroup$

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}, 
         Method -> {{Automatic, "LocalAdaptive"}[[-1]], 
           "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]] 

enter image description here

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

enter image description here

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

enter image description here

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

Hope it helps!

$\endgroup$
5
  • $\begingroup$ 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! $\endgroup$ Commented Jul 10, 2023 at 10:26
  • $\begingroup$ 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 $\endgroup$ Commented Jul 11, 2023 at 8:00
  • $\begingroup$ @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. $\endgroup$ Commented Jul 11, 2023 at 13:13
  • $\begingroup$ @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 $\endgroup$ Commented Jul 11, 2023 at 13:24
  • $\begingroup$ @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. $\endgroup$ Commented Jul 11, 2023 at 13:30
6
$\begingroup$

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}, 
          Method -> {{Automatic, "LocalAdaptive"}[[-1]], 
            "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]}]    

Figure 1 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}, 
          Method -> {{Automatic, "LocalAdaptive"}[[-1]], 
            "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]}] 

Figure 2

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]}] 

Figure 3

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}]]]

Figure 4

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}, 
        Method -> {{Automatic, "LocalAdaptive"}[[-1]], 
          "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. Figure 5

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}, 
         Method -> {{Automatic, "LocalAdaptive"}[[-1]], 
           "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}, 
        Method -> {{Automatic, "LocalAdaptive"}[[-1]], 
          "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[%]

Figure 6

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}, 
         Method -> {{Automatic, "LocalAdaptive"}[[-1]], 
           "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]]  
 

Figure 7

$\endgroup$
10
  • 1
    $\begingroup$ Yours second code dosen't work for r=2 ? I have no Plot ?I tried with MMA 13.3. $\endgroup$ Commented Jul 8, 2023 at 8:39
  • $\begingroup$ @MariuszIwaniuk Thank you. Don't forget to add parameter L=20; before plot. Anyway post has been updated using With[]. $\endgroup$ Commented Jul 8, 2023 at 9:52
  • $\begingroup$ Ok works fine now. Thanks $\endgroup$ Commented Jul 8, 2023 at 9:54
  • $\begingroup$ @AlexTrounev Nice answer!!! Does iterative solution also jump for case r==0 ? $\endgroup$ Commented Jul 8, 2023 at 10:36
  • $\begingroup$ @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. $\endgroup$ Commented Jul 8, 2023 at 12:09
5
$\begingroup$

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}]

enter image description here

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}}]

enter image description here

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}]

enter image description here

$\endgroup$
9
  • 1
    $\begingroup$ Your first Differentiate covers only the case r==0. This case has a general solution 1/(1+Exp[x])! $\endgroup$ Commented Jul 10, 2023 at 13:11
  • $\begingroup$ Note that for the $r = 2$ case the equation simply (?) becomes $a'''(x) = -2 e^{a(-x)}$. $\endgroup$ Commented Jul 11, 2023 at 14:45
  • $\begingroup$ @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 $\endgroup$ Commented Jul 11, 2023 at 16:41
  • 2
    $\begingroup$ @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 $\endgroup$ Commented Jul 12, 2023 at 8:04
  • 1
    $\begingroup$ @yarchik This is nice solution (+1). $\endgroup$ Commented Jul 14, 2023 at 7:16

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.