1
$\begingroup$

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 enter image description here

enter image description here

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}]
$\endgroup$
6
  • 3
    $\begingroup$ NDSolve cannot solve problems involving DiracDelta! You have to approximate for example as a limit. $\endgroup$ May 5, 2022 at 9:47
  • 1
    $\begingroup$ 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. $\endgroup$
    – Artes
    May 5, 2022 at 9:49
  • 3
    $\begingroup$ @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}]). $\endgroup$
    – user64494
    May 5, 2022 at 10:13
  • $\begingroup$ 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 $\endgroup$ May 7, 2022 at 4:31
  • $\begingroup$ 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. $\endgroup$
    – user64494
    May 7, 2022 at 6:21

1 Answer 1

3
$\begingroup$

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 = 
 NDSolve`FiniteDifferenceDerivative[Derivative[2], grid, 
  "DifferenceOrder" -> 2]; fddf3 = 
 NDSolve`FiniteDifferenceDerivative[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]

Figure 1

$\endgroup$
9
  • $\begingroup$ 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? $\endgroup$ May 8, 2022 at 6:22
  • $\begingroup$ 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. $\endgroup$ May 8, 2022 at 6:27
  • $\begingroup$ 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. $\endgroup$ May 8, 2022 at 6:34
  • $\begingroup$ How "a" is fixing s=0?@AlexTrounev $\endgroup$ May 8, 2022 at 6:36
  • $\begingroup$ Please check that grid[[1]]+a=0 therefore s starts from 0 in definition hR[s] . $\endgroup$ May 8, 2022 at 6:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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