NDSolve Integro-differential equation

Question

I have hard time to solve numerically the following integro-differential equation:

ξ0 = 39;
λ0 = 20;
max = 500;
B = 0.1;
NDSolve[
    {
    A''[x] - 1/(2 λ0^2 ξ0) Integrate[A[x1] Exp[-((x - x1)/ξ0)], {x1, 0, max}] == 0, 
    A'[0] == B, A[max] == 0
    },
    A,
    {x, 0, max}
]

once I run Mathematica I get the errors:

NDSolve::idelay: Initial history needs to be specified for all variables for delay-differential equations.

NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.`.

Is there anybody that can help me?

Thank you very much,

Mattia

Answer 1

Seems, an analytical solution is possible.

ξ0 = 39;
λ0 = 20;
max = 500;
B = 1/10;

integrand = E^(1/39 (-x + x1)) A[x1];

eq = -(Integrate[integrand, {x1, 0, 500}]/31200) + 
        A''[x]

Indefinite integration over x of A''[x] yields A'[x] and integration inside the x1-integral with integration constant r (I don't show all intermediate results here)

A'[x] == 1/31200 Integrate[Integrate[integrand, x] + r, {x1, 0, 500}]

Separate integration of r, the other part is 39*A''[x]

Edit: Correction of sign error

A'[x] == 1/31200 Integrate[r, {x1, 0, 500}] - 39 A''[x]

(*   Derivative[1][A][x] == (5 r)/312 - 39 (A^′′)[x]   *)

Since you know A'[0], you get

Derivative[1][A][0] == (5 r)/312 - 39 (A^′′)[0] == 1/10

Second integratio over x yield A[x]

A[x] == 1/31200 Integrate[
    Integrate[(r - 39 E^(-(x/39) + x1/39) A[x1]), x] + s, {x1, 0, 500}]

The s and r term is 5/312 (s + r x) plus 1521*A''[x]

1/31200 Integrate[s + r x, {x1, 0, 500}]

At x==500 you have

A[500] == 5/312 (500 r + s) + 1521 (A^′′)[500] == 0

Solve for r and s

sol1 = First@
     Solve[{(5 r)/312 - 39 A''[0] == 1/10, 
       5/312 (500 r + s) + 1521 A''[500] == 0}, {r, s}]

The differential equation is now eq2, which can be solved with DSolve

eq2 = A[x] == 5/312 (s + r x) + 1521 A''[x] /. sol1 // Simplify

Solve deq

dsol1 = First@
         DSolve[eq2 /. {A''[0] -> ass0, A''[500] -> ass500}, A, x]

(*   {A -> Function[{x}, 
     1/10 (-500 - 195000 ass0 - 15210 ass500 + x + 390 ass0 x) + 
     E^(x/39) C[1] + E^(-x/39) C[2]]}   *)

To eliminate C1 and C2 solve with boundary conditions

sol2 = First@
      Solve[{(A[500] /. dsol1) == 0, (A'[0] /. dsol1) == 1/10}, {C[1], 
       C[2]}]

now you have still a dependance of ass0 and ass500

A''[x] /. dsol1 /. sol2 // Simplify

(*   (E^(-x/39) (ass0 (E^(1000/39) - E^(2 x/39)) + 
      ass500 (E^(500/39) + E^((2 (250 + x))/39))))/(1 + E^(1000/39))   *)

Solve for ass0 and ass500 with the found function A

sol3 = First@
     Solve[{(A''[500] /. dsol1 /. sol2) == 
  ass500, (A''[0] /. dsol1 /. sol2) == ass0}, {ass500, ass0}] // 
   Simplify

(*   {ass0 -> ass500 E^(500/39)}   *)

Get remainig ass500 by comparing both sides of the equation

ls = A''[x] /. dsol1 /. sol2 /. sol3 // Simplify

rs = Integrate[integrand /. dsol1 /. sol2 /. sol3, {x1, 0, 500}]/31200

sol4 = First@Solve[ls == rs, ass500] // Simplify

(*   {ass500 -> -((539 - 39 E^(500/39))/(
      15210 + 382000 E^(500/39) - 15210 E^(1000/39)))}   *)

The desired function is then

A[x] /. dsol1 /. sol2 /. sol3 /. sol4 // Simplify[#, x > 0] &

(*   (E^(-x/39) (819819 E^(500/39) - 59319 E^(1000/39) + 
      E^((500 + x)/39) (8648819 - 17179 x) - 
      1521 E^(x/39) (39 + x)))/(10 (-1521 - 38200 E^(500/39) + 
      1521 E^(1000/39)))   *)

Test all conditions

A[500] /. dsol1 /. sol2 /. sol3 /. sol4 // Simplify[#, x > 0] &

(*   0   *)

A'[0] /. dsol1 /. sol2 /. sol3 /. sol4 // Simplify[#, x > 0] &

(*   1/10   *)

eq /. dsol1 /. sol2 /. sol3 /. sol4 // Simplify[#, x > 0] &

(*   0   *)

LogPlot[Evaluate[{-A[x], A[x]} /. dsol1 /. sol2 /. sol3 /. sol4 // 
    Simplify[#, x > 0] &], {x, 0, 500}, PlotStyle -> {Red, Blue}]

Plot[Evaluate[
  A[x] /. dsol1 /. sol2 /. sol3 /. sol4 // Simplify[#, x > 0] &], {x, 
  0, 500}, PlotRange -> All]

Answer 2

The $x$-dependent part of your integrand can be removed from the integral, leaving:

ode = A''[x] - Exp[-x/ξ0]/(2 λ0^2 ξ0) Integrate[A[t] Exp[t/ξ0], {t, 0, max}] == 0;

If we let:

b'[t] == A[t] Exp[t/ξ0]
b[0] == 0

then b[max] is equal to the integral. Let int be the value of the integral for the solution to your differential equation. Then, we expect the solution $A(x)$ to satisfy:

A''[x] - Exp[-x/ξ0]/(2 λ0^2 ξ0) int == 0

So, we are looking for the value of int where the above equation is satisfied, and b[max] == int. We can use ParametricNDSolveValue and FindRoot to do this:

pf = ParametricNDSolveValue[
    {
    A''[x] - Exp[-x/ξ0]/(2λ0^2 ξ0) int == 0, A'[0]==B, A[500]==0,
    b'[x] == A[x] Exp[x/ξ0], b[0]==0
    },
    {A,b[max]},
    {x,0,max},
    int
];

integral = i /. FindRoot[Indexed[pf[i], 2] == i, {i, 1}]

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

-80.0025

(I used Indexed instead of Part since Part will issue error messages for symbolic i)

Now that we know what the value of the integral is, we can determine A:

sol = pf[integral][[1]];

Visualization:

Plot[sol[t], {t, 0, max}, PlotRange->All]

Finally, here is a plot of the error:

Plot[sol''[x] - Exp[-x/ξ0]/(2λ0^2 ξ0) integral, {x, 0, 500}, PlotRange->All]

My results agree with @Akku's.

Answer 3

The following is a FDM approach whose result is agree with Akku14 and Carl Woll's. I've used pdetoae for the generation of difference equation.

ξ0 = 39;
λ0 = 20;
max = 500;
B = 1/10;

SetAttributes[int, Listable];
eq = A''[x] - 1/(2 λ0^2 ξ0) int[x] == 0;
kernel[x_, x1_] = A[x1] Exp[-((x - x1)/ξ0)];
bc = {A'[0] == B, A[max] == 0};

points = 25;
difforder = 4;
domain = {0, max};

{nodes, weights} = Most[NIntegrate`GaussRuleData[points, MachinePrecision]];
midgrid = Rescale[nodes, {0, 1}, domain];


intrule = int@x_ :> -Subtract @@ domain weights.Map[kernel[x, #] &, midgrid];

grid = Flatten[{domain // First, midgrid, domain // Last}];
(* Definition of pdetoae isn't included in this post,
   please find it in the link above. *)
ptoafunc = pdetoae[A[x], grid, difforder];
del = #[[2 ;; -2]] &;
ae = del@ptoafunc[eq] /. intrule;
aebc = ptoafunc@bc;
(*initialguess[x_]=-10;
sollst=FindRoot[{ae,aebc},Table[{A@x,initialguess@x},{x,grid}]][[All,-1]];*)
sollst = Solve[{ae, aebc} // Flatten, A /@ grid][[1, All, -1]];
sol = Interpolation[{grid, sollst}\[Transpose]];

Plot[sol@x, {x, 0, max}, PlotRange -> All]

---

Update

If you feel the usage of del confused, the followings are 2 alternatives that don't require one to remove redudant equations:

fullae = ptoafunc[eq] /. intrule;

(* Approach 1 *)
lSSolve[obj_List, constr___, x_, opt : OptionsPattern[FindMinimum]] := 
 FindMinimum[{1/2 obj^2 // Total, constr}, x, opt]
lSSolve[obj_, rest__] := lSSolve[{obj}, rest]

sollst = lSSolve[Subtract @@@ Flatten[{fullae, aebc}], A /@ grid][[2, All, -1]];

(* Approach 2 *)
{blst, mat} = CoefficientArrays[Flatten@{fullae, aebc}, A /@ grid];
sollst = LeastSquares[N@mat, -blst];
sol = Interpolation[{grid, sollst}\[Transpose]];

Plot[sol@x, {x, 0, max}, PlotRange -> All]

If you want to learn more about lSSolve, check this post.