3
$\begingroup$

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.3188535755423771`Exp[-0.07679009872276484`(x-x1)])-4.529403364443408`*^-17/(1+1.3148412969171435`Exp[-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]?

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

1 Answer 1

Reset to default
3
$\begingroup$

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]

Figure 1

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 Figure 2 We can combine all plots in one using Show Figure 3

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

Figure 4

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

Figure 5

Improve this answer
$\endgroup$
13
  • $\begingroup$ Thank you very much for your answer, Alex. I will review it and let you know. $\endgroup$
    – Asatur Khurshudyan
    May 15, 2022 at 10:43
  • 1
    $\begingroup$ How did you choose the values for k0 and M0? Both affect the solution pretty much. $\endgroup$
    – Asatur Khurshudyan
    May 15, 2022 at 23:40
  • 1
    $\begingroup$ 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/… $\endgroup$
    – Alex Trounev
    May 16, 2022 at 2:08
  • $\begingroup$ I can recommend to add some constraint like f>0 at 0<x<l. $\endgroup$
    – Alex Trounev
    May 16, 2022 at 3:30
  • 1
    $\begingroup$ 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. $\endgroup$
    – Alex Trounev
    May 16, 2022 at 4:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.