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.

$\endgroup$
4
  • 1
    $\begingroup$ Can you add the definition of g[tau] to your code? $\endgroup$ Sep 16 at 21:00
  • $\begingroup$ You're right, just added it! $\endgroup$ Sep 17 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$ Sep 18 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$ Sep 20 at 8:52
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

$\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$ Sep 24 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$ Sep 24 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$ Sep 25 at 15:57
  • $\begingroup$ @JoseEnriqueAroca May be it could better to combine curves in two plots? See update to my answer. $\endgroup$ Sep 26 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

$\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$ Sep 17 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$ Sep 20 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 at 9:30

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.