7
$\begingroup$

I am trying to solve the following coupled differential equations: $\dot{c}_1(t)=-\int_0^tf(t-s)c_1(s)+g(t-s)c_2(s)ds$, $\dot{c}_2(t)=-\int_0^tg(t-s)c_1(s)+f(t-s)c_2(s)ds$, where $f(t-s)=\frac{\sqrt{2}\pi^{\frac{3}{2}}g_{ab}^2\rho_0r_0^2e^{i\omega(t-s)}}{\hbar^2\left(r_0^2+\frac{i\hbar(t-s)}{m_b}\right)^{\frac{3}{2}}}$ and $g(t-s)=\frac{\sqrt{2}e^{i\omega(t-s)-\frac{\Delta r_c^2m_b}{2m_br_0^2+2i(t-s)\hbar}}g_{ab}^2m_b^{\frac{3}{2}}\pi^{\frac{3}{2}}r_0^2\rho_0\left(-\Delta r_c^2m_b+m_br_0^2+i(t-s)\hbar\right)}{\hbar^2\left(m_br_0^2+i\hbar(t-s)\right)^{\frac{5}{2}}}$.

I developed a code that uses the finite difference method to solve an equation in the form $\dot{c}_1(t)=-\int_0^tf(t-s)c_1(s)ds$, but I am unsure how to extend this to solve the coupled differential equations. My attempt to do so is as follows:

Clear[gab, \[Rho]0, \[HBar], \[Omega]0, mb, \[Tau], k, r0, t, s, drc]
drc = 0;
\[HBar] = 1;
\[Omega]1 = 10;
\[Omega]0 = 10;
gab = 1;
mb = 1;
ma = 1;
\[Rho]0 = 1;
drc = 0;
r0 = Sqrt[(\[HBar])/(ma*\[Omega]0)];
f[\[Tau]_] := (
  Sqrt[2] E^(I (\[Tau]) \[Omega]0) gab^2 \[Pi]^(3/2)
    r0^2 \[Rho]0)/(\[HBar]^2 (r0^2 + (I (\[Tau]) \[HBar])/mb)^(
   3/2)) ;
g[\[Tau]_] := (Sqrt[2] E^(
     I \[Omega]0 \[Tau] - (drc^2 mb)/(
      2 mb r0^2 + 2 I (\[Tau]) \[HBar])) gab^2 mb^(3/2) \[Pi]^(3/2)
      r0^2 \[Rho]0 (-drc^2 mb + mb r0^2 + 
       I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb r0^2 + 
       I (\[Tau]) \[HBar])^(5/2));
dt = 0.1;
nsubint = 1000;
ds = dt/nsubint;
Ttab = Table[T, {T, 0, 10, dt}];
Stab = Table[s, {s, 0, dt - ds, ds}];
(*Ti[T_]:=Position[Ttab,T][[1]]*)
Clear[Ti];
Clear[c1, c2];
c1[0] = 1;
c2[0] = 0;
Clear[cTtab1, cTab2];
(*c1[1]=c1[0]-dt(c1[0]*f[dt]+c2[0]*g[dt];*)
\
(*c2[1]=c2[0]-dt(c1[0]*g[dt]+c2[0]*f[dt];*)
Do[
 corrSum1[n] = 
  Sum[c1[nn - 1]*
     Sum[f[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn , 
       1}], {nn, 1, n}] + 
   Sum[c2[nn - 1]*
     Sum[g[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn , 
       1}], {nn, 1, n}];
 corrSum2[n] = 
  Sum[c1[nn - 1]*
     Sum[g[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn , 
       1}], {nn, 1, n}] + 
   Sum[c2[nn - 1]*
     Sum[f[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn , 
       1}], {nn, 1, n}];
 c1[n] = c1[n - 1] - dt*corrSum1[n]; 
 c2[n] = c2[n - 1] - dt*corrSum2[n], {n, 0, 100}]
cTtab1 = Table[{n*dt, c1[n]}, {n, 0, 100}]
cTtab2 = Table[{n*dt, c2[n]}, {n, 0, 100}]
FDiff1 = ListPlot[Abs[cTtab1], PlotStyle -> Orange, 
  PlotLegends -> {"Finite Difference c1"}]
FDiff2 = ListPlot[Abs[cTtab2], PlotStyle -> Orange, 
  PlotLegends -> {"Finite Difference c2"}]

The absolute value of c1 and c2 should go between 0 and 1.

Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ Can you add the definition of g[tau] to your code? $\endgroup$
    – Michael Seifert
    Sep 16, 2021 at 21:00
  • $\begingroup$ You're right, just added it! $\endgroup$
    – Jose Enrique Aroca
    Sep 17, 2021 at 12:44
  • $\begingroup$ This is a system of integrodifferential equations. Since f, g are complex we have 4 real equations as it shown by Cesario. $\endgroup$
    – Alex Trounev
    Sep 18, 2021 at 8:16
  • $\begingroup$ @JoseEnriqueAroca Your code is quite well, but it is not so precise as, for instance, iterative method (see my answer). $\endgroup$
    – Alex Trounev
    Sep 20, 2021 at 8:52

2 Answers 2

Reset to default
7
$\begingroup$

We can use the original code after small correction as follows

Clear[gab, \[Rho]0, \[HBar], \[Omega]0, mb, \[Tau], k, r0, t, s, drc]
drc = 0;
\[HBar] = 1;
\[Omega]1 = 10;
\[Omega]0 = 10;
gab = 1;
mb = 1;
ma = 1;
\[Rho]0 = 1;
drc = 0;
r0 = Sqrt[(\[HBar])/(ma*\[Omega]0)];
f[\[Tau]_] := (Sqrt[
      2] E^(I (\[Tau]) \[Omega]0) gab^2 \[Pi]^(3/
        2) r0^2 \[Rho]0)/(\[HBar]^2 (r0^2 + (I (\[Tau]) \[HBar])/
         mb)^(3/2));
g[\[Tau]_] := (Sqrt[
      2] E^(I \[Omega]0 \[Tau] - (drc^2 mb)/(2 mb r0^2 + 
           2 I (\[Tau]) \[HBar])) gab^2 mb^(3/2) \[Pi]^(3/
        2) r0^2 \[Rho]0 (-drc^2 mb + mb r0^2 + 
       I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb r0^2 + 
        I (\[Tau]) \[HBar])^(5/2));
dt = 0.01;
nsubint = 10;
ds = dt/nsubint;
Ttab = Table[T, {T, 0, 10, dt}];
Stab = Table[s, {s, 0, dt - ds, ds}];
(*Ti[T_]:=Position[Ttab,T][[1]]*)
Clear[c1, c2];
c1[0] = 1; c1[-1] = 1;
c2[0] = 0; c2[-1] = 0;
Clear[cTtab1, cTtab2];
Do[corrSum1[n] = 
  Sum[c1[nn - 1]*
     Sum[f[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn, 
       1}], {nn, 1, n}] + 
   Sum[c2[nn - 1]*
     Sum[g[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn, 
       1}], {nn, 1, n}];
 corrSum2[n] = 
  Sum[c1[nn - 1]*
     Sum[g[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn, 
       1}], {nn, 1, n}] + 
   Sum[c2[nn - 1]*
     Sum[f[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn, 
       1}], {nn, 1, n}];
 c1[n] = c1[n - 1] - dt*corrSum1[n];
 c2[n] = c2[n - 1] - dt*corrSum2[n], {n, 0, 100}]

cTtab1 = Table[{n*dt, Abs[c1[n]]}, {n, 0, 100}];
cTtab2 = Table[{n*dt, Abs[c2[n]]}, {n, 0, 100}];

FDiff1 = ListPlot[cTtab1, PlotStyle -> Orange, 
  PlotLegends -> {"Finite Difference c1"}]
FDiff2 = ListPlot[cTtab2, PlotStyle -> Orange, 
  PlotLegends -> {"Finite Difference c2"}]

Figure 1

The second method is iteration with substitution s->t/2 (s+1), and therefore Volterra equations are reduced to Fredholm equations as follows

X[0][t_] := 1; 
Y[0][t_] := 0; ds = 
 1/100; nmax = 15; Do[{X[i], Y[i]} = 
   NDSolveValue[{x'[
       t] == -t ds/2 Sum[
        f[t - t/2 (s + 1)] X[i - 1][t/2 (s + 1)] + 
         g[t - t/2 (s + 1)] Y[i - 1][t/2 (s + 1)], {s, -1, 1, ds}], 
     y'[t] == -t ds/2 Sum[
        g[t - t/2 (s + 1)] X[i - 1][t/2 (s + 1)] + 
         f[t - t/2 (s + 1)] Y[i - 1][t/2 (s + 1)], {s, -1, 1, ds}], 
     x[0] == 1, y[0] == 0}, {x, y}, {t, 0, 1}];, {i, 1, nmax}]

Visualization of FDM and iterative solution

{Show[Plot[Evaluate[Abs[X[nmax][t]]], {t, 0, 1}, 
   FrameLabel -> {t, "c1"}, Frame -> True], FDiff1], 
 Show[Plot[Evaluate[Abs[Y[nmax][t]]], {t, 0, 1}, 
   FrameLabel -> {t, "c2"}, Frame -> True], FDiff2]}

Figure 2

We also can compare the iterative method with code @Cesario. After small modification we have

Clear["Global`*"]
n = 100;
tmax = 1;
dT = tmax/n; drc = 0;
\[HBar] = 1;
\[Omega]1 = 10;
\[Omega]0 = 10;
gab = 1;
mb = 1;
ma = 1;
\[Rho]0 = 1;
drc = 0;
r0 = Sqrt[(\[HBar])/(ma*\[Omega]0)];
f[\[Tau]_] := (Sqrt[
      2] E^(I (\[Tau]) \[Omega]0) gab^2 \[Pi]^(3/
        2) r0^2 \[Rho]0)/(\[HBar]^2 (r0^2 + (I (\[Tau]) \[HBar])/
         mb)^(3/2));
g[\[Tau]_] := (Sqrt[
      2] E^(I \[Omega]0 \[Tau] - (drc^2 mb)/(2 mb r0^2 + 
           2 I (\[Tau]) \[HBar])) gab^2 mb^(3/2) \[Pi]^(3/
        2) r0^2 \[Rho]0 (-drc^2 mb + mb r0^2 + 
       I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb r0^2 + 
        I (\[Tau]) \[HBar])^(5/2));
fr[k_] := Re[N[f[k dT]]]
fi[k_] := Im[N[f[k dT]]]
gr[k_] := Re[N[g[k dT]]]
gi[k_] := Im[N[g[k dT]]]

vars = Flatten[Table[{c1r[k], c1i[k], c2r[k], c2i[k]}, {k, 0, n}]];
equ1r[k_] := (c1r[k] - c1r[k - 1])/dT + 
  dT Sum[(-c1i[j] fi[1 - j + k] + c1r[j] fr[1 - j + k] - 
      c2i[j] gi[1 - j + k] + c2r[j] gr[1 - j + k]), {j, 1, k}]
equ1i[k_] := (c1i[k] - c1i[k - 1])/dT + 
  dT Sum[(c1r[j] fi[1 - j + k] + c1i[j] fr[1 - j + k] + 
      c2r[j] gi[1 - j + k] + c2i[j] gr[1 - j + k]), {j, 1, k}]
equ2r[k_] := (c2r[k] - c2r[k - 1])/dT + 
  dT Sum[(-c2i[j] fi[1 - j + k] + c2r[j] fr[1 - j + k] - 
      c1i[j] gi[1 - j + k] + c1r[j] gr[1 - j + k]), {j, 1, k}]
equ2i[k_] := (c2i[k] - c2i[k - 1])/dT + 
  dT Sum[(c2r[j] fi[1 - j + k] + c2i[j] fr[1 - j + k] + 
      c1r[j] gi[1 - j + k] + c1i[j] gr[1 - j + k]), {j, 1, k}]
equs = Flatten[
   Table[{equ1r[k] == 0, equ1i[k] == 0, equ2r[k] == 0, 
     equ2i[k] == 0}, {k, 1, n}]];
equstot = 
  Join[equs, {c1r[0] - 1 == 0, c1i[0] == 0, c2r[0] == 0, c2i[0] == 0}];
sol = Solve[equstot, vars][[1]];



p1 = ListPlot[Table[{k dT, Norm[{c1r[k], c1i[k]}] /. sol}, {k, 0, n}],
   PlotStyle -> Red]

p2 = ListPlot[Table[{k dT, Norm[{c2r[k], c2i[k]}] /. sol}, {k, 0, n}],
   PlotStyle -> Red]

Comparison of iterative method and Cesario's method (X, Y are taken from the code above)

{Show[Plot[Evaluate[Abs[X[nmax][t]]], {t, 0, 1}, 
   FrameLabel -> {t, "c1"}, Frame -> True], p1], 
 Show[Plot[Evaluate[Abs[Y[nmax][t]]], {t, 0, 1}, 
   FrameLabel -> {t, "c2"}, Frame -> True], p2]}

Figure 3

Note, that all 3 methods give us qualitatively same results, but they are differ in detail.
For parametric research we can use Module, for example,

plot[om_] := Module[{\[Omega]0 = om, f, g}, drc = 0;
  \[HBar] = 1;
  \[Omega]1 = 10;
  gab = 1;
  mb = 1;
  ma = 1;
  \[Rho]0 = 1;
  drc = 0;
  r0 = Sqrt[(\[HBar])/(ma*\[Omega]0)];
  f[\[Tau]_] := (Sqrt[
       2] E^(I (\[Tau]) \[Omega]0) gab^2 \[Pi]^(3/
         2) r0^2 \[Rho]0)/(\[HBar]^2 (r0^2 + (I (\[Tau]) \[HBar])/
          mb)^(3/2));
  g[\[Tau]_] := (Sqrt[
       2] E^(I \[Omega]0 \[Tau] - (drc^2 mb)/(2 mb r0^2 + 
            2 I (\[Tau]) \[HBar])) gab^2 mb^(3/2) \[Pi]^(3/
         2) r0^2 \[Rho]0 (-drc^2 mb + mb r0^2 + 
        I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb r0^2 + 
         I (\[Tau]) \[HBar])^(5/2));
  dt = 0.01;
  nsubint = 10;
  ds = dt/nsubint;
  Ttab = Table[T, {T, 0, 10, dt}];
  Stab = Table[s, {s, 0, dt - ds, ds}];
  Clear[c1, c2];
  c1[0] = 1; c1[-1] = 1;
  c2[0] = 0; c2[-1] = 0;
  Clear[cTtab1, cTtab2];
  Do[corrSum1[n] = 
    Sum[c1[nn - 1]*
       Sum[f[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn, 
         1}], {nn, 1, n}] + 
     Sum[c2[nn - 1]*
       Sum[
        g[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn, 1}], {nn,
        1, n}];
   corrSum2[n] = 
    Sum[c1[nn - 1]*
       Sum[g[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn, 
         1}], {nn, 1, n}] + 
     Sum[c2[nn - 1]*
       Sum[f[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn, 
         1}], {nn, 1, n}];
   c1[n] = c1[n - 1] - dt*corrSum1[n];
   c2[n] = c2[n - 1] - dt*corrSum2[n], {n, 0, 100}];
  cTtab1 = Table[{n*dt, Abs[c1[n]]}, {n, 0, 100}];
  cTtab2 = Table[{n*dt, Abs[c2[n]]}, {n, 0, 100}]; {cTtab1, cTtab2}]

With this code we can plot solution as function of $\omega 0$ as follows

pl = Table[
  ListPlot[plot[om], PlotRange -> All, PlotLegends -> {"c1", "c2"}, 
   PlotLabel -> Row[{"\[Omega]0 = " , om}], 
   AxesLabel -> {"t", ""}], {om, {5, 10, 15}}]

Figure 4

We can also combine curves in two plots

pl1 = ListPlot[Table[plot[om][[1]], {om, {5, 10, 15}}], 
  PlotRange -> All, AxesLabel -> {"t", "c1"}, 
  PlotLegends -> {5, 10, 15}]

pl2 = ListPlot[Table[plot[om][[2]], {om, {5, 10, 15}}], 
  PlotRange -> All, AxesLabel -> {"t", "c2"}, 
  PlotLegends -> {5, 10, 15}]

Figure 5

Improve this answer
$\endgroup$
4
  • $\begingroup$ Hi, thank you very much for your detailed answer. I have a small question. Say I want to use the original code to do some plots with changing parameters (e.g. different Omega_1 values) and plot the results on the same axis. What would be the most efficient way to do this? Can I create an array of omega_2 values or something similar? $\endgroup$
    – Jose Enrique Aroca
    Sep 24, 2021 at 11:15
  • $\begingroup$ @JoseEnriqueAroca Your code not depends on $\omega 1$, but depends on $\omega 0$. We can use your code as Module to plot solution vs different $\omega 0$ - see update to my answer. $\endgroup$
    – Alex Trounev
    Sep 24, 2021 at 16:44
  • $\begingroup$ Hi yeah, sorry for the typo, I meant 𝜔0. Is there any way to show the plots for the three 𝜔0 values on the same axis (i.e show the 6 lines on the same axis) $\endgroup$
    – Jose Enrique Aroca
    Sep 25, 2021 at 15:57
  • $\begingroup$ @JoseEnriqueAroca May be it could better to combine curves in two plots? See update to my answer. $\endgroup$
    – Alex Trounev
    Sep 26, 2021 at 3:31
2
$\begingroup$

To have a coarse view of $c_1,c_2$ we can proceed as follows

Clear["Global`*"]
n = 100;
tmax = 5;
dT = tmax/n;
rc = 1;
\[HBar] = 1;
\[Omega] = 10;
gab = 1;
mb = 1;
ma = 1;
\[Rho]0 = 1;
r0 = 1;
\[CapitalDelta] = 1;
f[k_] := (E^(I k dT \[Omega]) gab^2 \[Pi]^(3/2) r0^2\[Rho]0)/(Sqrt[2] \[HBar]^2 (r0^2 + (I k dT \[HBar])/mb)^(3/2));
g[k_] := (Sqrt[2] E^(I k dT \[Omega] - \[CapitalDelta] rc^2 mb/(2 mb r0^2 + 2 I k dT \[HBar])) gab^2 mb^(3/2) \[Pi]^(3/2) r0^2 \[Rho]0) (-\[CapitalDelta] rc^2 mb + mb r0^2 + I k dT \[HBar])/( \[HBar]^2 (mb r0^2 + (I k dT \[HBar]))^(5/2));

fr[k_] := Re[N[f[k]]]
fi[k_] := Im[N[f[k]]]
gr[k_] := Re[N[g[k]]]
gi[k_] := Im[N[g[k]]]

vars = Flatten[Table[{c1r[k], c1i[k], c2r[k], c2i[k]}, {k, 0, n}]];
equ1r[k_] := (c1r[k] - c1r[k - 1])/dT + dT Sum[(- c1i[j] fi[1 - j + k] + c1r[j] fr[1 - j + k] - c2i[j] gi[1 - j + k] + c2r[j] gr[1 - j + k]), {j, 1, k}]
equ1i[k_] := (c1i[k] - c1i[k - 1])/dT + dT Sum[( c1r[j] fi[1 - j + k] + c1i[j] fr[1 - j + k] + c2r[j] gi[1 - j + k] + c2i[j] gr[1 - j + k]), {j, 1, k}]
equ2r[k_] := (c2r[k] - c2r[k - 1])/dT + dT Sum[(- c2i[j] fi[1 - j + k] + c2r[j] fr[1 - j + k] - c1i[j] gi[1 - j + k] + c1r[j] gr[1 - j + k]), {j, 1, k}]
equ2i[k_] := (c2i[k] - c2i[k - 1])/dT + dT Sum[( c2r[j] fi[1 - j + k] + c2i[j] fr[1 - j + k] + c1r[j] gi[1 - j + k] + c1i[j] gr[1 - j + k]), {j, 1, k}]
equs = Flatten[Table[{equ1r[k], equ1i[k], equ2r[k], equ2i[k]}, {k, 1, n}]];
equstot = Join[equs, {c1r[0] - 1/2, c1i[0], c2r[0], c2i[0] - 1/2}];
sol = Solve[equstot == 0, vars][[1]];

(* ListPlot[Table[{k dT, c1r[k]}, {k, 0, n}] /. sol]
ListPlot[Table[{k dT, c1i[k]}, {k, 0, n}] /. sol] *)
ListPlot[Table[{c1r[k], c1i[k]}, {k, 0, n}] /. sol]
(* ListPlot[Table[{k dT, c2r[k]}, {k, 0, n}] /. sol]
ListPlot[Table[{k dT, c2i[k]}, {k, 0, n}] /. sol] *)
ListPlot[Table[{c2r[k], c2i[k]}, {k, 0, n}] /. sol]

enter image description here

enter image description here

Improve this answer
$\endgroup$
3
  • $\begingroup$ Hi, thank you very much for your answer. I think that this would work fine. However, I have been asked to explicitly use the finite difference method shown above. I have edited the code in the question with my attempt to include c1 and c2 in the same Do Loop, but for some reason, I get answers for c1 and c2 bigger than 1 (they should never exceed this value). You can take a look if you like. $\endgroup$
    – Jose Enrique Aroca
    Sep 17, 2021 at 12:46
  • $\begingroup$ Nice code (+1), but you have typo in definition of fr[k_] := Re[N[f[k]]]; fi[k_] := Im[N[f[k]]]; gr[k_] := Re[N[g[k]]]; gi[k_] := Im[N[g[k]]]. It should be fr[k_] := Re[N[f[k dT]]];... Please, have a look Figure 3 in my answer with comparison your method and iterative method. $\endgroup$
    – Alex Trounev
    Sep 20, 2021 at 8:50
  • $\begingroup$ Thanks for your appreciation (+1). In my case, the k*dT transformation is included into the function definition. $\endgroup$
    – Cesareo
    Sep 20, 2021 at 9:30

Your Answer

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

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