# Numerical solution of an advance-delay ODE with multiple advances and delays

I am trying to solve a differential equation of the form

a[x_] := x
b[x_] := x^2
c1[x_] := x
x1 = 1;
c2[x_] := x
x2 = 10;
c3[x_] := x
x3 = 5;
c4[x_] := x
x4 = 15;
eq = D[a[x] f'[x] + b[x] f[x], x] + (c1[x + x1] f[x + x1] - c1[x] f[x]) + (c2[x + x2] f[x + x2] - c2[x] f[x]) + (c3[x - x3] f[x - x3] - c3[x] f[x]) + (c4[x - x4] f[x - x4] - c4[x] f[x]) == 0;


with the boundary conditions at "infinity"

l = 100;
bc1 = f'[l] == 0;
bc2 = f[l] == 0;


Real coefficients are way to difficult and the number of advance-delay terms is larger, so this is just an example problem.

All coefficients a, b, c1, ..., c4 are zero everywhere outside the integration region [0, l]. Moreover, we can assume that f is also 0 outside [0, l].

I know that NDSolve handles delay and advance differential equations pretty well. But I can't figure out how to apply it for solving this system. Any help is appreciated.

Edit

I should have mentioned that the unknown function f[x] satisfies the integral constraint

Integrate[Sqrt[x] f[x], {x, 0, l}] == 1


Update

I am updating the OP with some corrections to parameters and specific expressions for the system coefficients.

First, I have changed xmin and therefore the definition of l, so that the integration stops exactly at x = 1.

k0 =3;
M0 =8;
xmin = .025;
xmax = 100;
l = xmax;
x1 = 11.65;


Limiting the consideration by only c1 term, I have obtained the coefficients of the real system as follows:

dAfFun[x_]:=Piecewise[{{1.4476805109594064 10^-23+1.1735563659975857 10^-24  x^(1.4),xmin<=x<=xmax}},0]
AfFun[x_]:=Piecewise[{{4.17425813460343 10^-24+4.4423429301477025 10^-23 x,xmin<=x<=xmax}},0]
Rx1Fun[x_]:=Piecewise[{{4.537932520459608 10^-17/(1+1.3188535755423771Exp[-0.07679009872276484(x-x1)])-4.529403364443408*^-17/(1+1.3148412969171435Exp[-0.07683234639209759(x-x1)]),x1<=x<=xmax}},0]


Thence,

a[x_] := dAfFun[x l]
b[x_] := AfFun[x l]
c1[x_] := Rx1Fun[x l]


As a result, NSolve and FindRoot find negative roots (the solution looks reasonable though), while NMinimize with Method -> "RandomSearch" runs for unreasonably long time.

The actual solution pretty much looks like PDF[ChiDistribution[1], x] plus effects from the c1 term visible after x1. So, maybe instead of EulerE, we need to consider PDF[ChiDistribution[k], x]?

• What have you tried so far? Commented May 11, 2022 at 1:35
• I have tried to solve it with a single delay and single advance separately. But I am not able to combine them. Commented May 11, 2022 at 4:07
• With boundary conditions bc1 = f'[l] == 0; bc2 = f[l] == 0 solution is $f(x)=0$. Commented May 11, 2022 at 17:04
• @AsaturKhurshudyan With additional constraint you added there is nonzero solution - see my answer. Commented May 14, 2022 at 8:10
• @AsaturKhurshudyan Please, see Update 2 to my answer. Commented May 17, 2022 at 11:51

This problem can be solved by colocation method with using Euler wavelets. First, we map solution on unit interval $$0\le x\le 1$$. Then we define function f, first (f1) and second (f2) derivative as follows

UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 7; xmin = 0; xmax = 100; l = xmax - xmin; nn =
Total[With[{k = k0, M = M0},
Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]]; dx =
1/(nn); xl = Table[xmin + s*dx, {s, 0, nn}]; xcol =
Table[(xl[[s - 1]] + xl[[s]])/2, {s, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Intx1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Intx2 = Integrate[Intx1, t1];
intx1[y_] := Intx1 /. t1 -> y; intx2[y_] := Intx2 /. t1 -> y;
Psi[y_] := Psijk /. t1 -> y;

a[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
b[x_] := l^2 Piecewise[{{x^2, 0 <= x <= 1}, {0, True}}];
c1[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x1 = 1/l;
c2[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x2 = 10/l;
c3[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x3 = 5/l;
c4[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x4 = 15/l;
(*eq=D[a[x] f'[x]+b[x] f[x],x]+(c1[x+x1] f[x+x1]-c1[x] \
f[x])+(c2[x+x2] f[x+x2]-c2[x] f[x])+(c3[x-x3] f[x-x3]-c3[x] \
f[x])+(c4[x-x4] f[x-x4]-c4[x] f[x])\[Equal]0;bc1=f'[l]\[Equal]0;
bc2=f[l]\[Equal]0;*)

var = Array[v, {nn}]; f[x_] := var . intx2[x] + x v1 + v0;
f1[x_] := var . intx1[x] + v1; f2[x_] := var . Psi[x];


Finally we compute equation in collocation points and solve system of algebraic equations with FindRoot

eq = Table[
a[x] f2[x]/l^2 + f1[x]/l + b[x] f1[x]/l +
2 l x f[x] + (c1[x + x1] f[x + x1] -
c1[x] f[x]) + (c2[x + x2] f[x + x2] -
c2[x] f[x]) + (c3[x - x3] f[x - x3] -
c3[x] f[x]) + (c4[x - x4] f[x - x4] - c4[x] f[x]) == 0, {x,
xcol}];

eqs = Join[
Drop[eq, -1], {l^1.5 Sum[
Sqrt[xcol[[i]]] f[xcol[[i]]] (xl[[i + 1]] - xl[[i]]), {i,
nn}] == 1, f[1] == 0, f1[1] == 0}]; varM = Join[var, {v0, v1}];

sol = FindRoot[eqs, Table[{varM[[i]], 1/10}, {i, Length[varM]}]]


Visualization f[x]

Plot[f[x/l] /. sol, {x, 0, l}, Frame -> True, PlotRange -> All]


Update 1. We also can use NMinimize with bound f[x]>0 at 0<x<l as follows (note we define exact integral intf[x] instead of sum)

UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 8; xmin = 0; xmax = 100; l = xmax - xmin; nn =
Total[With[{k = k0, M = M0},
Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]]; dx =
1/(nn); xl = Table[xmin + s*dx, {s, 0, nn}]; xcol =
Table[(xl[[s - 1]] + xl[[s]])/2, {s, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Intx1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Intx2 = Integrate[Intx1, t1]; Intx3 =
Integrate[Intx2 Sqrt[t1], t1, Assumptions -> t1 > 0];
intx1[y_] := Intx1 /. t1 -> y; intx2[y_] := Intx2 /. t1 -> y;
intx3[y_] := Intx3 /. t1 -> y; Psi[y_] := Psijk /. t1 -> y;

a[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
b[x_] := l^2 Piecewise[{{x^2, 0 <= x <= 1}, {0, True}}];
c1[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x1 = 1/l;
c2[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x2 = 10/l;
c3[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x3 = 5/l;
c4[x_] := l Piecewise[{{x, 0 <= x <= 1}, {0, True}}];
x4 = 15/l;
(*eq=D[a[x] f'[x]+b[x] f[x],x]+(c1[x+x1] f[x+x1]-c1[x] \
f[x])+(c2[x+x2] f[x+x2]-c2[x] f[x])+(c3[x-x3] f[x-x3]-c3[x] \
f[x])+(c4[x-x4] f[x-x4]-c4[x] f[x])\[Equal]0;bc1=f'[l]\[Equal]0;
bc2=f[l]\[Equal]0;*)

var = Array[v, {nn}];
intf[x_] := var . intx3[x] + 2/5 v1 x^(5/2) + 2/3 v0 x^(3/2);
f[x_] := var . intx2[x] + x v1 + v0; f1[x_] := var . intx1[x] + v1;
f2[x_] := var . Psi[x];

eq = Table[
a[x] f2[x]/l^2 + f1[x]/l + b[x] f1[x]/l +
2 l x f[x] + (c1[x + x1] f[x + x1] -
c1[x] f[x]) + (c2[x + x2] f[x + x2] -
c2[x] f[x]) + (c3[x - x3] f[x - x3] -
c3[x] f[x]) + (c4[x - x4] f[x - x4] - c4[x] f[x]), {x, xcol}];

bc =
Join[{l^(3/2) intf[1] == 1, f[1] == 0, f1[1] == 0},
Table[f[x] > 0, {x, xcol}]]; varM = Join[var, {v0, v1}];
solm =
NMinimize[{eq . eq, bc}, varM]


Visualization of solutions for nn=10, 16, 24, 28, 32, 56 We can combine all plots in one using Show

Update 2. We can use wavelets generated by $$e^{-m t}$$ to solve second problem as follows

UE[m_, t_] := Exp[-m t]
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 7; nn =
Total[With[{k = k0, M = M0},
Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]]; dx =
1/(nn); xl = Table[s*dx, {s, 0, nn}]; xcol =
Table[(xl[[s - 1]] + xl[[s]])/2, {s, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Intx1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Intx2 = Integrate[Intx1 Sqrt[t1], t1, Assumptions -> t1 > 0];
intx1[y_] := Intx1 /. t1 -> y; intx2[y_] := Intx2 /. t1 -> y;
Psi[y_] := Psijk /. t1 -> y;

xmin = 0; xmax = 100; l = xmax - xmin; x0 = 1; x1 = x0/l;
dAfFun[x_] :=
Piecewise[{{1.4476805109594064 10^-23 +
1.1735563659975857 10^-24 x^(1.4), xmin <= x <= xmax}}, 0]
AfFun[x_] :=
Piecewise[{{4.17425813460343 10^-24 + 4.4423429301477025 10^-23 x,
xmin <= x <= xmax}}, 0]
Rx1Fun[x_] :=
Piecewise[{{4.537932520459608 10^-17/(1 +
1.3188535755423771 Exp[-0.07679009872276484 (x - x0)]) -
4.529403364443408*^-17/(1 +
1.3148412969171435 Exp[-0.07683234639209759 (x - x0)]),
x0 <= x <= xmax}}, 0]

a[x_] := dAfFun[x l] 10^20
b[x_] := AfFun[x l] 10^20
c1[x_] := Rx1Fun[x l] 10^20

var1 = Array[v, {nn}]; var2 = Array[u, {nn}];
intf[x_] := var1 . intx2[x] + 2/3 v0 x^(3/2);
f[x_] := var1 . intx1[x] + v0; f1[x_] := var1 . Psi[x];
g1[x_] := var2 . Psi[x]; g[x_] := var2 . intx1[x] + u0;

(*eq={a[x] f'[x]+b[x] f[x]==g[x],D[g[x],x]+(c1[x+x1] f[x+x1]-c1[x] \
f[x])==0};bc1=f'[l]\[Equal]0;
bc2=f[l]\[Equal]0;*)
eq =
Join[Table[g1[x]/l + (c1[x + x1] f[x + x1] - c1[x] f[x]), {x, xcol}],
Table[-g[x] + a[x] f1[x]/l + b[x] f[x], {x, xcol}]]; varM =
Join[var1, var2, {v0, u0}];
eqs = Join[
Table[eq[[i]] == 0, {i, 2, Length[eq] - 1}], {l^(3/2) intf[1] == 1,
f1[0] == 0, f[1] == 0, g[1] == 0}]; sol =
FindRoot[eqs, Table[{varM[[i]], 1/10}, {i, Length[varM]}]];


Visualization numerical solution and PDF[ChiDistribution[1], x]

{Plot[PDF[ChiDistribution[1], x], {x, 0, 10}, PlotRange -> All,
Frame -> True, PlotLabel -> "PDF[ChiDistribution[1],x]"],
Plot[Evaluate[{f[x/l] /. sol}], {x, 0, 20}, Frame -> True,
PlotRange -> All, PlotLabel -> nn],
Plot[Evaluate[{f[x/l] /. sol}], {x, 0, 100}, Frame -> True,
PlotRange -> All, PlotLabel -> nn]}


Note, that we can solve linear system of equations with LinearSolve using

{vec, mat} = CoefficientArrays[eqs, varM];

soll = LinearSolve[mat, -vec];

rule = Table[varM[[i]] -> soll[[i]], {i, Length[varM]}];


Visualization two numerical solutions

Plot[Evaluate[{f[x/l] /. rule, f[x/l] /. sol}], {x, 0, l},
Frame -> True, PlotRange -> All, PlotLabel -> nn,
PlotLegends -> {"LinearSolve", "FindRoot"},
PlotStyle -> {Blue, Dashed}]


• Thank you very much for your answer, Alex. I will review it and let you know. Commented May 15, 2022 at 10:43
• How did you choose the values for k0 and M0? Both affect the solution pretty much. Commented May 15, 2022 at 23:40
• It is why we need some theorem about method convergence for this equation. For pantograph type delay differential equations this method working well also with Haar wavelets - see, for instance, sciencedirect.com/science/article/pii/… Commented May 16, 2022 at 2:08
• I can recommend to add some constraint like f>0 at 0<x<l. Commented May 16, 2022 at 3:30
• This is starting point for iterations. I have computed several runs with nn=10,16, 24, 28, 32, 56 and then combine in one plot. It looks like we need to restrict f[0] as well, for example, f[0]==0 or f[0]==1. See update 1 to my answer. Commented May 16, 2022 at 4:43