# Issues with convergence of an iteration of numerical integrals computed with NIntegrate

This is my first question on Stack-Exchange. I am not an expert in numerical integration and I am having issues of convergence with some numerical integration. I am trying to solve numerically a nonlinear integral equation of which I know an exact solution. I want to compare my numerical results with the exact one (to check the correctness of the code). Unfortunately, after two iterations the numerical result starts departing from the exact behavior.

The integral equation to solve is:

$$\qquad y(t)=\int_{-\infty}^{\infty}dx G(t-x) [y(x)^2-1]y(x)$$,

where

$$\qquad G(t-x) = (1/2) (t - x)\, {\rm sign}(t - x).$$

An analytic solution is given by

$$\qquad y(x)=\tanh(x/\sqrt{2})$$

The code I am using is the following:

t1 = -10;
t2 = 10;
n = 100;
m = 5;
t = Range[t1, t2, (t2 - t1)/(n)];
Do[
y[x_, 0] = Tanh[x/Sqrt[2]];
list[j] =
{t,
NIntegrate[
(1/2)(t - x) Sign[t-x](y[x, j - 1]^2 - 1) y[x, j-1], {x, -10, 10},
AccuracyGoal -> 13, WorkingPrecision -> 13]} // Transpose;
sol[j] = Interpolation[list[j], Method -> "Spline"];
y[x_, j] = sol[j][x],
{j, m}]


Basically, I am trivially starting from the exact solution for j = 0, and check the convergence at higher order in the iteration (I am doing this to verify whether the code is correct). But from j = 3 the numerical solution starts departing from the exact one; so I would like to find a way to correct my code and make the iteration convergent for any higher j, as it should be.

I hope I have clearly stated my problem.

• I don't know much about Green's Function, but doesn't it require a linear differential operator? Commented Oct 22, 2020 at 16:44
• Hi Josh. Basically, given the equation y''(t)=F[y(t)], we can always formally write the solution as an integral of G(t-t') * F[y(t')], where G is the green function defined as the solution of the differential equation with a Dirac delta source. But I am not able to get numerically a convergent solution for that integral, whose exact solution is known. Any help would be very much appreciated. Thanks. Commented Oct 22, 2020 at 17:23
• It looks like you try to develop a new method of solution of nonlinear Fredholm integral equations. Have you any theorem for this solution convergence? Commented Oct 22, 2020 at 17:56
• @user92 If I now understand you right you're looking for an iterative solution of the integralequation? If yes, starting the iteration with the known exact solution evaluates to True , that means the iteration will stop. Commented Oct 23, 2020 at 7:25
• @user92 Yes, find a "good" starting function and start iteration with NestList Commented Oct 23, 2020 at 7:40

Not a complete answer, but an attempt to derive the correct GreenFunction and the equivalent integralequation using greens method.

The basic idea of @user92 is the (mis)use of greenfunction for a nonlinear ode.

Greenfunction for the linear ode y''[x]==f[x] is given by

g[x_, \[Xi]_] := (x - \[Xi]) HeavisideTheta[x - \[Xi]]



!!!Attention, @user92 gave another possible wrong formula!!!

A possible solution of the linear ode y''[x]==f[x] can be expressed as y[x]=Integrate[(x - \[Xi]) f[\[Xi]], {\[Xi], -Infinity, x}]

D[Integrate[(x - \[Xi]) f[\[Xi]], {\[Xi], -Infinity, x}], {x, 2}]
(*f[x]*)


If this intergal holds for arbitrary f[x] it's valid to substitude the right hand side of the nonlinear ode y''[x]== (y[x]^2-1)y[x]

This gives the integral equation

Integrate[(x - \[Xi]) (y[\[Xi]]^2 - 1) y[\[Xi]],{\[Xi], -Infinity,x}] == y[x]

The given exact solution y[x]=Tanh[x/Sqrt[2]] isn't causal and fullfills the modified integralequation Integrate[(x - \[Xi]) (y[\[Xi]]^2 - 1) y[\[Xi]],{\[Xi], -Infinity,x}]-1 == y[x]

Y = -1 + Integrate[(x - \[Xi]) (y[\[Xi]]^2 - 1) y[\[Xi]], {\[Xi], -Infinity, x}] /.y -> ( Tanh[#/Sqrt[2]] &) // Simplify
Plot[{Y , Tanh[x/Sqrt[2]]}, {x, -5, 5},PlotStyle -> {Thickness[.01], Dashed}]


An iterative solution of the integralequation might be found using NestList

• Since in your nonlinear Volterra integral equation $f=-1$, it is why numerical solution converses to $\tanh(x/\sqrt{2})$ not to $y=0$. Commented Oct 23, 2020 at 11:32
• @AlexTrounev Yes agreed, but it's not so easy to find a starting function for the iteration! Commented Oct 23, 2020 at 11:46
• The starting function it is also a reason why people prefer collocation method. Commented Oct 23, 2020 at 12:24
• @AlexTrounev But the method needs correct(!) integralequation, developped from the given nonlinear ode, to work and give the right solution! Commented Oct 23, 2020 at 13:13
• @UlrichNeumann: Many thanks for your answer. I think that the way I wrote the question was a bit confusing, so I edited it trying to make it less confusing. With Green function I just meant the kernel of the integral equation, but I agree that it was a bit confusing. I now hope that my problem is clearly stated. Commented Oct 23, 2020 at 16:01

It is a common problem for iterative method of solution of integral equations that numerical error decreases in the few first iteration and then increases. In this sense we prefer standard collocation method like described here. For example, we have nonlinear Fredholm integral equation $$u(x)=f(x)+\int_0^1G(t,x)(u(t)^2-1)u(t)dt$$
with $$G(t,x)=1/2 (t-x) sign(t-x)$$, and with exact solution $$u(x)=\tanh (x/\sqrt{2})$$. Integrated exact solution we get that

f[t_]:=-(1/4) (-Sqrt[2] t + Sqrt[2] Sech[1/Sqrt[2]]^2 -
Sqrt[2] t Sech[1/Sqrt[2]]^2 - 2 Tanh[1/Sqrt[2]]);


For numerical solution we use Bernoulli wavelets and GaussianQuadratureWeights[]. This code produces numerical solution with absolute error $$2.26\times 10^{-5}$$ with 32 collocation points only:

Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"]; \
Get["NumericalDifferentialEquationAnalysis"];
ue[t_] := Tanh[t/Sqrt[2]];
f[t_] := -(1/
4) (-Sqrt[2] t + Sqrt[2] Sech[1/Sqrt[2]]^2 -
Sqrt[2] t Sech[1/Sqrt[2]]^2 - 2 Tanh[1/Sqrt[2]]);
G[x_, z_] := (1/2) (x - z) Sign[x - z];
n = 4;
M = Sum[1, {j, 0, n, 1}, {i, 0, 2^j - 1, 1}] + 1;
dx = 1/M; A = 0; xl = Table[A + l*dx, {l, 0, M}]; xcol =
Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, M + 1}];
psi1[x_] := Piecewise[{{BernoulliB[2, x], 0 <= x < 1}, {0, True}}];
psi2[x_] :=  Piecewise[{{BernoulliB[1, x], 0 <= x < 1}, {0, True}}];
psi1jk[x_, j_, k_] := psi1[j*x - k];
psi2jk[x_, j_, k_] := psi2[j*x - k];
psijk[x_, j_, k_] := (psi1jk[x, j, k] + psi2jk[x, j, k]);

np = 2 M; g = GaussianQuadratureWeights[np, 0, 1];
points = g[[All, 1]]; weights = g[[All, 2]];
GuassInt[ff_, z_] :=
Sum[(ff /. z -> points[[i]])*weights[[i]], {i, np}];
u[t_] := Sum[
a[j, k]*psijk[t, 2^j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}] + a0 ;
int[x_] := GuassInt[G[x, z] (u[z]^2 - 1) u[z], z];
eq = Table[-u[xcol[[i]]] + f[xcol[[i]]] + int[xcol[[i]]] == 0, {i,
Length[xcol]}];
varM = Join[{a0},
Flatten[Table[a[j, k], {j, 0, n, 1}, {k, 0, 2^j - 1, 1}]]];
sol = FindRoot[eq, Table[{varM[[i]], 1/10}, {i, Length[varM]}]];
unum = Table[ {xcol[[i]], Evaluate[u[xcol[[i]]] /. sol]}, {i,
Length[xcol]}];

(*Absolute error*)  du = Table[{x, Abs[ue[x] - Evaluate[u[x] /. sol]]}, {x, xcol}]


Visualization

{Show[Plot[ue[x], {x, 0, 1},
PlotLegends ->
Placed[LineLegend[{"Exact"}, LabelStyle -> {Black, 15}],
Scaled[{0.2, 0.8}]], AspectRatio -> 1,
LabelStyle -> Directive[{FontSize -> 15}, Black],
AxesLabel -> {"x", "u"}, PlotStyle -> Blue],
ListPlot[unum, PlotRange -> All, PlotStyle -> Red,
PlotLegends ->
Placed[PointLegend[{"Numeric"}, LabelStyle -> {Black, 15}],
Scaled[{0.2, 0.9}]]]],
ListPlot[du, Filling -> Axis, PlotRange -> All,
PlotLabel ->
Row[{"Absolute error. Number of collocation points =", M}]]}


• If the integralequation holds especially for f[x]==0, the expression u[x] - NIntegrate[G[x, z] (u[z]^2 - 1) u[z], {z, 0, 1}] /. u -> (Tanh[#/Sqrt[2]] &) should vanish for all x` I think. But it doesn't. Any idea? Commented Oct 23, 2020 at 6:51
• @Alex Trounev: many thanks for your answer, but I still have further questions. I need to solve numerically the integral with f[x]=0 and with the interval from -\infinity to +\infinity, because then I want to apply the same method to another integral equation with a different kernel and of which I don't know any exact solution. So, do you think that your code can be adapted to the case f[x]=0 and interval from -infty to +infty (of course by taking something like -a to +a, where a is sufficiently large). Commented Oct 23, 2020 at 6:59
• @UlrichNeumann Hello! I think that the that difference (with f[x]=0) will vanish only when the integral is taken from -infinity to +infinity. I mean that the hyperbolic tangent is an exact solution of the integral with f[x]=0 only when we take the integral from -infinity to +infinity. Commented Oct 23, 2020 at 7:06
• @user92 Thanks, I checked it, you're right! That means, you should change the integration range in your question too? Commented Oct 23, 2020 at 7:11
• @UlrichNeumann: You are right! I am sorry, I didn't specify it....I will change it now. I had put from -10 to 10 meaning some sufficiently large interval. Thanks for the notice. Commented Oct 23, 2020 at 7:17