# Solving integro-differential Equation with Dirac-Delta Term For Dynamical Analysis

I am trying to solve Equation.44 from this paper https://labsites.rochester.edu/agrawal/wp-content/uploads/2019/08/paper_2019_03.pdf But I am unable to tackle the Integrand term which includes Dirac-Delta term and an additional response function as follows

Here I attach my code in MATHEMATICA. Is there any process so that I can check the normalization throughout the program? As the particle number should be 3, I took U_0=1.5. So that after normalization I will get N=3.

a = 20;
Cvalue = 0.5;
qvalue = 100;
Nvalue = Sqrt[Cvalue];
Fvalue[ζ_] :=
Cos[π*qvalue*ζ]^2 + Cvalue^2 Sin[π*qvalue*ζ]^2;
para1 = 12.2*10^-15;
para2 = 32*10^-15;
fR = 0.18;
δ = 0.02;
hR[k_] := (para1^2 + para2^2)/(para1*para2^2)
Exp[-k/para2] Sin[k/para1];
R[k_] := (1 - fR) DiracDelta[k] + fR*hR[k];
eqs = {I D[u[ζ, τ], ζ] +
1/2 D[u[ζ, τ], {τ, 2}] -
I*δ*D[u[ζ, τ], {τ, 3}] +
Nvalue^2/
Fvalue[ζ] u[ζ, τ] Integrate[
R[s] Abs[u[ζ, τ - s]]^2, {s, 0, ∞}] == 0,
u[0, τ] == 1.5 Sech[τ], u[ζ, -a] == 0,
u[ζ, a] == 0};
sol = NDSolveValue[eqs, u, {τ, -a, a}, {ζ, 0, 1},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 20, "MaxPoints" -> 81,
"DifferenceOrder" -> "Pseudospectral"}}];
DensityPlot[
Abs[sol[ζ, τ]]^2, {τ, -a, a}, {ζ, 0, 1},
PlotPoints -> 100, ColorFunction -> "Rainbow"]


So, I omitted the Dirac-Delta term and still failed here is the code

a = 20; \[Zeta]max = 1; step = 10^-6; d\[Tau] = (a - step)/6 - step;
n = 5;
int[0][\[Zeta]_, \[Tau]_] := 0;
Cvalue = 0.5;
qvalue = 100;
Nvalue = Sqrt[Cvalue];
Fvalue[\[Zeta]_] :=
Cos[\[Pi]*qvalue*\[Zeta]]^2 + Cvalue^2 Sin[\[Pi]*qvalue*\[Zeta]]^2;
\[Delta] = 0.02;
fR = 0.18;
Do[Itaretion[i] =
NDSolveValue[{I D[u[\[Zeta], \[Tau]], \[Zeta]] +
1/2 D[u[\[Zeta], \[Tau]], {\[Tau], 2}] -
I \[Delta] D[u[\[Zeta], \[Tau]], {\[Tau], 3}] +
Nvalue^2/Fvalue[\[Zeta]] u[\[Zeta], \[Tau]]*fR*
int[i - 1][\[Zeta], \[Tau]] == 0,
u[\[Zeta], -L] == u[\[Zeta], L],
u[0, \[Tau]] == 1.5*Sech[\[Tau]]},
u, {\[Tau], -L, L}, {\[Zeta], 0, \[Zeta]max},
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 40}}}];
int[i] =
Interpolation[
NIntegrate[
Abs[U[i][\[Zeta], \[Tau] - s]]^2*Exp[-s/10^-15] Sin[s/10^-15], {s,
0, \[Infinity]}], {\[Tau], -a, a, d\[Tau]}, {\[Zeta],
0, \[Zeta]max, .2*\[Zeta]max}];, {i, 1, n}]

DensityPlot[
Itaretion[n][\[Zeta], \[Tau]], {\[Tau], -a, a}, {\[Zeta],
0, \[Zeta]max}]

• NDSolve cannot solve problems involving DiracDelta! You have to approximate for example as a limit. May 5, 2022 at 9:47
• You should write down your equations in the question, explain what they describe and clarify your code. If you don't do it I will down vote your question. May 5, 2022 at 9:49
• @UlrichNeumann: Your statement "NDSolve cannot solve problems involving DiracDelta" does not correspond to reality (see the result of sol = NDSolve[{f'[x] == f[x] + DiracDelta[x - 1], f[0] == 1}, f, {x, 0, 2}]; Plot[Evaluate[f[x]/.sol],{x,0,2}]). May 5, 2022 at 10:13
• I have tried solving the equation omitting only the Dirac-Delta term in the integration and allowing the h_R(t) term. But still I am unable to do it.@UlrichNeumann@user64494 May 7, 2022 at 4:31
• I looked in the linked article. It should be noticed that the notation $\delta(x-x')$ appears (for the first time) in (16) without any description and/or definition . The statement "that the Kernel in Eq. (11) is reduced to the form" before (16) is unclear to me. May 7, 2022 at 6:21

This problem can be solved with using method of lines as here. First, we note that $$\int_{s_1}^{s_2} \delta(s)f(t-s)=f(t)$$ for any interval with $$s1<0,s_2>0$$. But if we take $$s_1=0, s_2\rightarrow\infty$$ as in a paper, then we have not well defined result

Integrate[DiracDelta[s] f[t - s], {s, 0, Infinity}]

Out[]= -f[t] (-1 + HeavisideTheta[0])


Therefore the answer depends on how we interpreter HeavisideTheta[0]. In numerical model it looks like HeavisideTheta[0]=RandomReal[], or may be 0.5, or may be 0. In our computation we put HeavisideTheta[0]=0. Second question is about function $$h_R(t)$$ normalization. In our computation we use for normalization $$T_0=57 fs$$. Finally we have

a = 20; h = 1/10; grid = Range[-a, a, h]; nn = Length[grid];
varu[t_] := Table[u[i][t], {i, nn}]; var =
Table[u[i], {i, 2, nn - 1}]; fddf2 =
NDSolveFiniteDifferenceDerivative[Derivative[2], grid,
"DifferenceOrder" -> 2]; fddf3 =
NDSolveFiniteDifferenceDerivative[Derivative[3], grid,
"DifferenceOrder" -> 2]; mat2 =
fddf2@"DifferentiationMatrix"; mat3 =
fddf3@"DifferentiationMatrix"; du2 = mat2 . varu[t]; du3 =
mat3 . varu[t];
Cvalue = .5;
qvalue = 100;
Nvalue = 3;
Fvalue[t_] := Cos[\[Pi]*qvalue*t]^2 + Cvalue^2 Sin[\[Pi]*qvalue*t]^2;
T0 = 57(*fs*);
para1 = 12.2/T0;
para2 = 32/T0;
fR = 0.18;
\[Delta] = 0.02;
hR[k_] := (para1^2 + para2^2)/(para1*para2^2) Exp[-k/para2] Sin[
k/para1];
R[k_] :=(*(1-fR) DiracDelta[k]+fR*)hR[k];
eqs = Table[
I D[u[i][t], t] + 1/2 du2[[i]] - I*\[Delta]*du3[[i]] +
Nvalue^2/
Fvalue[t] ((1 - fR) Abs[u[i][t]]^2 +
fR h Sum[
If[i - s < 1, 0, hR[grid[[s]] + a] Abs[u[i - s][t]]^2], {s,
1, nn}]) u[i][t] == 0, {i, 3, nn - 1}];
ic = Table[u[i][0] == Sech[grid[[i]]], {i, 3, nn - 1}]; u[1][t_] := 0;
u[nn][t_] := 0; u[2][t_] := 0

sol = NDSolve[{eqs, ic}, var, {t, 0, 1}, MaxSteps -> 10^6,
Method -> "ExplicitRungeKutta"];


Visualization

lst = Table[{{t, grid[[i]]}, Abs[u[i][t]] /. sol[[1]]}, {t, 0,
1, .01}, {i, nn}]; U = Interpolation[Flatten[lst, 1]];
Plot3D[U[t, x], {t, 0, 1}, {x, -a, a}, PlotRange -> All,
PlotPoints -> 50, Mesh -> None, ColorFunction -> Hue,
PlotTheme -> "Marketing", AxesLabel -> {\[Zeta], \[Tau]}]

DensityPlot[U[t, x], {x, -a, a}, {t, 0, 1}, PlotRange -> All,
PlotPoints -> 50, ColorFunction -> Hue,
AxesLabel -> {\[Tau], \[Zeta]}, PlotLegends -> Automatic]


• Thanks, Sir. I have a couple of questions. 1. The reason for Multiplying h before summing in eqs. and 2. Addition of "a" in If command. Is there any reason to use the Explicite Rungekutta method? May 8, 2022 at 6:22
• This is approximation of $\int_0^{\infty} h_R(s)|u(\zeta ,\tau-s)|^2ds$ with $ds=h$. Adding a allow us to fix s=0. No reason it just example of Method. May 8, 2022 at 6:27
• Option  MaxSteps -> 10^6 used for the case with {t,0,100}. In a case {t,0,1} we can remove all options in NDSolve. May 8, 2022 at 6:34
• How "a" is fixing s=0?@AlexTrounev May 8, 2022 at 6:36
• Please check that grid[[1]]+a=0 therefore s starts from 0 in definition hR[s] . May 8, 2022 at 6:42