2
$\begingroup$

I am trying to solve the following coupled differential equations: $\dot{c}_{1}(t)=-\int_0^t \alpha_{1}(t,s)c_{1}(s)+\beta_{1}(t,s)c_{2}(s)\,ds-\frac{i}{\hbar}\lambda c_2(t)e^{+i(\omega_1-\omega_2) t}$ $\dot{c}_{2}(t)=-\int_0^t \alpha_{2}(t,s)c_{2}(s)+\beta_{2}(t,s)c_{1}(s)\,ds-\frac{i}{\hbar}\lambda c_1(t)e^{-i(\omega_1-\omega_2) t}$

I developed a code that uses the finite difference method to solve the equations $\dot{c}_{1}(t)=-\int_0^t \alpha_{1}(t,s)c_{1}(s)+\beta_{1}(t,s)c_{2}(s)\,ds$
$\dot{c}_{2}(t)=-\int_0^t \alpha_{2}(t,s)c_{2}(s)+\beta_{2}(t,s)c_{1}(s)\,ds$,
where $\alpha_1(t,s)=\frac{\sqrt{2}\pi^{\frac{3}{2}}g_{ab}^2\rho_0r_{01}^2e^{i\omega_{1}(t-s)}}{\hbar^2\left(r_{01}^2+\frac{i\hbar(t-s)}{m_b}\right)^{\frac{3}{2}}}$,
$\alpha_2(t,s)=\frac{\sqrt{2}\pi^{\frac{3}{2}}g_{ab}^2\rho_0r_{02}^2e^{i\omega_{2}(t-s)}}{\hbar^2\left(r_{02}^2+\frac{i\hbar(t-s)}{m_b}\right)^{\frac{3}{2}}}$,
$\beta_1(t,s)=\frac{4\pi^{\frac{3}{2}}g_{ab}^2\rho_0r_{01}r_{02}e^{i\omega_{1}(t-s)}m_b^{\frac{3}{2}}e^{i(\omega_1-\omega_2)s}e^{-\frac{\Delta r_c^2m_b}{m_b\left(r_{01}^2+r_{02}^2\right)+2i(t-s)\hbar}}}{{\hbar^2\left(m_b\left(r_{01}^2+r_{02}^2\right)+2i(t-s)\hbar\right)^{\frac{5}{2}}}}$
and $\beta_2(t,s)=\frac{4\pi^{\frac{3}{2}}g_{ab}^2\rho_0r_{01}r_{02}e^{i\omega_{2}(t-s)}m_b^{\frac{3}{2}}e^{-i(\omega_1-\omega_2)s}e^{-\frac{\Delta r_c^2m_b}{m_b\left(r_{01}^2+r_{02}^2\right)+2i(t-s)\hbar}}}{{\hbar^2\left(m_b\left(r_{01}^2+r_{02}^2\right)+2i(t-s)\hbar\right)^{\frac{5}{2}}}}$
but I am unsure how to extend this to include the $-\frac{i}{\hbar}\lambda c_{12}(t)e^{+-i(\omega_1-\omega_2) t}$ factors (lambda is just a constant). My code to solve the first part of the equations is as follows.

Clear[gab, \[Rho]0, \[HBar], mb, r01, r02, drc, \[Omega]1, \[Omega]2, \
\[Omega]]
drc = 2*r01;
\[HBar] = 1;
\[Omega]1 = 50;
\[Omega]2 = 55;
\[Omega] = (\[Omega]1*\[Omega]2)/(\[Omega]1 + \[Omega]2);
gab = 1;
mb = 1;
ma = 1;
\[Rho]0 = 1;
r01 = Sqrt[(\[HBar])/(ma*\[Omega]1)];
r02 = Sqrt[(\[HBar])/(ma*\[Omega]2)];
\[Alpha]1[\[Tau]_] := (
  Sqrt[2] E^(I (\[Tau]) \[Omega]1) gab^2 \[Pi]^(3/2)
    r01^2 \[Rho]0)/(\[HBar]^2 (r01^2 + (I (\[Tau]) \[HBar])/mb)^(
   3/2)) ;
\[Alpha]2[\[Tau]_] := (
  Sqrt[2] E^(I (\[Tau]) \[Omega]2) gab^2 \[Pi]^(3/2)
    r02^2 \[Rho]0)/(\[HBar]^2 (r02^2 + (I (\[Tau]) \[HBar])/mb)^(
   3/2)) ;
\[Beta]1[\[Tau]_] := (4 E^(
    I  \[Omega]1 \[Tau] - (drc^2 mb)/(
     mb (r01^2 + r02^2) + 2 I (\[Tau]) \[HBar])) gab^2 mb^(
    3/2) \[Pi]^(3/2)
     r01 r02 \[Rho]0 (-2 drc^2 mb + mb (r01^2 + r02^2) + 
      2 I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb (r01^2 + r02^2) + 
      2 I (\[Tau]) \[HBar])^(5/2))
\[Beta]2[\[Tau]_] := (4 E^(
    I \[Omega]2 \[Tau] - (drc^2 mb)/(
     mb (r01^2 + r02^2) + 2 I (\[Tau]) \[HBar])) gab^2 mb^(
    3/2) \[Pi]^(3/2)
     r01 r02 \[Rho]0 (-2 drc^2 mb + mb (r01^2 + r02^2) + 
      2 I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb (r01^2 + r02^2) + 
      2 I (\[Tau]) \[HBar])^(5/2))

dt = 0.015;
nsubint = 1000;
ds = dt/nsubint;
Clear[c1, c2];
c1[0] = 1;
c2[0] = 0;
dp = 30;
\[Lambda] = 10;

Clear[cTtab1, cTtab2];
Do[
 corrSum1[n] = 
  Sum[c1[nn - 1]*
     Sum[\[Alpha]1[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn ,
        1}], {nn, 1, n}] + 
   Sum[c2[nn - 1]*
     Sum[E^(I (\[Omega]1 - \[Omega]2) (m ds)) \[Beta]1[n dt - m ds]*
       ds, {m, nsubint (nn - 1), nsubint nn , 1}], {nn, 1, n}];
 corrSum2[n] = 
  Sum[c1[nn - 1]*
     Sum[E^(I (\[Omega]2 - \[Omega]1) (m ds)) \[Beta]2[n dt - m ds]*
       ds, {m, nsubint (nn - 1), nsubint nn , 1}], {nn, 1, n}] + 
   Sum[c2[nn - 1]*
     Sum[\[Alpha]2[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn ,
        1}], {nn, 1, n}];
 c1[n] = c1[n - 1] - dt*corrSum1[n]-dt*I*\[Lambda]*c2[n - 1];
 c2[n] = c2[n - 1] - dt*corrSum2[n]-dt*I*\[Lambda]*c1[n - 1], {n, 1, dp}]
cTtab1 = Table[{n*dt*\[Omega], Abs[c1[n]]^2}, {n, 0, dp}];
cTtab2 = Table[{n*dt*\[Omega], Abs[c2[n]]^2}, {n, 0, dp}];
Improve this question
$\endgroup$
11
  • $\begingroup$ It looks like we discussed this problem once on mathematica.stackexchange.com/questions/255697/… $\endgroup$
    – Alex Trounev
    Commented Nov 22, 2021 at 4:09
  • $\begingroup$ Hi, thank you very much for your comment! Indeed, we discussed a similar problem, but the introduction of a new function which is only dependent on t makes things more complicated. I thought that creating a new question would be more useful for new people to see it and add there ideas on how to solve it! $\endgroup$
    – Jose Enrique Aroca
    Commented Nov 22, 2021 at 8:59
  • $\begingroup$ Did you use code from my answer or this is a new code? $\endgroup$
    – Alex Trounev
    Commented Nov 22, 2021 at 16:13
  • $\begingroup$ This code is similar to what I posted previously. I finally didn't use the iterative solution you proposed. I managed to make a small adjustment to what I posted to make this working and now I have to include an extra factor, but as this is outside the integral, I am not sure how to. $\endgroup$
    – Jose Enrique Aroca
    Commented Nov 22, 2021 at 16:20
  • $\begingroup$ Did you test your code or you use it as it is? In my answer there are 3 methods are compared to test your code. $\endgroup$
    – Alex Trounev
    Commented Nov 23, 2021 at 3:36

1 Answer 1

Reset to default
3
$\begingroup$

We can solve this problem using FDM and implicit method of integration differential equations as follows:

Clear[gab, \[Rho]0, \[HBar], mb, r01, r02, drc, \[Omega]1, \[Omega]2, \
\[Omega]]
drc = 2*r01;
\[HBar] = 1; \[Lambda] = 10; k = I \[Lambda]/\[HBar];
\[Omega]1 = 50;
\[Omega]2 = 55; dom = \[Omega]2 - \[Omega]1;
\[Omega] = (\[Omega]1*\[Omega]2)/(\[Omega]1 + \[Omega]2);
gab = 1;
mb = 1;
ma = 1;
\[Rho]0 = 1;
r01 = Sqrt[(\[HBar])/(ma*\[Omega]1)];
r02 = Sqrt[(\[HBar])/(ma*\[Omega]2)];
f1[\[Tau]_] := (Sqrt[
      2] E^(I (\[Tau]) \[Omega]1) gab^2 \[Pi]^(3/
        2) r01^2 \[Rho]0)/(\[HBar]^2 (r01^2 + (I (\[Tau]) \[HBar])/
         mb)^(3/2));
f2[\[Tau]_] := (Sqrt[
      2] E^(I (\[Tau]) \[Omega]2) gab^2 \[Pi]^(3/
        2) r02^2 \[Rho]0)/(\[HBar]^2 (r02^2 + (I (\[Tau]) \[HBar])/
         mb)^(3/2));
g1[\[Tau]_] := (4 E^(I \[Omega]1 \[Tau] - (drc^2 mb)/(mb (r01^2 + 
              r02^2) + 2 I (\[Tau]) \[HBar])) gab^2 mb^(3/2) \[Pi]^(3/
        2) r01 r02 \[Rho]0 (-2 drc^2 mb + mb (r01^2 + r02^2) + 
       2 I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb (r01^2 + r02^2) + 
        2 I (\[Tau]) \[HBar])^(5/2));
g2[\[Tau]_] := (4 E^(I \[Omega]2 \[Tau] - (drc^2 mb)/(mb (r01^2 + 
              r02^2) + 2 I (\[Tau]) \[HBar])) gab^2 mb^(3/2) \[Pi]^(3/
        2) r01 r02 \[Rho]0 (-2 drc^2 mb + mb (r01^2 + r02^2) + 
       2 I (\[Tau]) \[HBar]))/(\[HBar]^2 (mb (r01^2 + r02^2) + 
        2 I (\[Tau]) \[HBar])^(5/2));

dt = 0.01;
nsubint = 10;
ds = dt/nsubint;
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[f1[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn, 
       1}], {nn, 1, n}] + 
   Sum[c2[nn - 1]*
     Sum[g1[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn, 
       1}], {nn, 1, n}];
 corrSum2[n] = 
  Sum[c1[nn - 1]*
     Sum[g2[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn, 
       1}], {nn, 1, n}] + 
   Sum[c2[nn - 1]*
     Sum[f2[n dt - m ds]*ds, {m, nsubint (nn - 1), nsubint nn, 
       1}], {nn, 1, n}];
 c1[n] = -((-c1[-1 + n] + 1/2 dt E^(I dom dt (-1 + n)) k c2[-1 + n] + 
    dt corrSum1[n] - 
    1/2 dt E^(I dom dt n)
      k (1/2 dt E^(-I dom dt (-1 + n)) k c1[-1 + n] - c2[-1 + n] + 
       dt corrSum2[n]))/(1 - (dt^2 k^2)/4)); 
 c2[n] = -(1/(-4 + dt^2 k^2)) E^(-I dom dt (-1 + n) - 
    I dom dt n) (-2 dt E^(I dom dt (-1 + n)) k c1[-1 + n] - 
     2 dt E^(I dom dt n) k c1[-1 + n] + 
     4 E^(I dom dt (-1 + n) + I dom dt n) c2[-1 + n] + 
     dt^2 E^(2 I dom dt (-1 + n)) k^2 c2[-1 + n] + 
     2 dt^2 E^(I dom dt (-1 + n)) k corrSum1[n] - 
     4 dt E^(I dom dt (-1 + n) + I dom dt n) corrSum2[n]);, {n, 0, 
  500}]

Visualization

cTtab1 = Table[{(n + 1)*dt, Abs[c1[n]]}, {n, 0, 500}];
cTtab2 = Table[{(n + 1)*dt, Abs[c2[n]]}, {n, 0, 500}];



pl1 = ListLinePlot[{cTtab1, cTtab2}, PlotRange -> All, 
  PlotLegends -> {"c1", "c2"}, AxesLabel -> {"t", ""}]

Figure 1

We can compare FDM result (points) with iterative method (solid lines) explained in my answer here

X[0][t_] := 1;
Y[0][t_] := 0; ds1 = 
 1/100; tmax = 1; nmax = 10; Do[{X[i], Y[i]} = 
   NDSolveValue[{x'[t]/
       tmax == -tmax t ds1/2 Sum[
         f1[tmax (t - t/2 (s + 1))] X[i - 1][t/2 (s + 1)] + 
          g1[tmax (t - t/2 (s + 1))] Y[i - 1][t/2 (s + 1)], {s, -1, 1,
           ds1}] - k y[t] Exp[I dom tmax t], 
     y'[t]/tmax == -tmax t ds1/2 Sum[
         g2[tmax (t - t/2 (s + 1))] X[i - 1][t/2 (s + 1)] + 
          f2[tmax (t - t/2 (s + 1))] Y[i - 1][t/2 (s + 1)], {s, -1, 1,
           ds1}] - k x[t] Exp[-I dom tmax t], x[0] == 1, 
     y[0] == 0}, {x, y}, {t, 0, 1}];, {i, 1, nmax}]

pl2 = Plot[Evaluate[{X[nmax][t], Y[nmax][t]} // Abs], {t, 0, 1}, 
  PlotLegends -> {"c1", "c2"}];
pl3 = ListPlot[{Take[cTtab1, 100], Take[cTtab2, 100]}, 
  PlotRange -> All];
Show[pl2, pl3]

Figure 2

Improve this answer
$\endgroup$
4
  • $\begingroup$ Hi, thank you for the great answer. I have some doubts. In the FDM, how have you accounted for the exp(i dom s) factor? And what is the motivation behind using c1[n]=-... and not just the positive? Also, is there a way to extend the iterative method beyond t=1.0? $\endgroup$
    – Jose Enrique Aroca
    Commented Nov 28, 2021 at 12:28
  • $\begingroup$ To get c1[n],c2[n] just evaluate Solve[c1[n] == c1[n - 1] - dt*corrSum1[n] - (c2[n - 1] Exp[I dom (n - 1) dt] + c2[n] Exp[I dom n dt])/2 dt k && c2[n] == c2[n - 1] - dt*corrSum2[n] - (c1[n - 1] Exp[-I dom (n - 1) dt] + c1[n] Exp[-I dom n dt])/2 dt k, {c1[n], c2[n]}]. $\endgroup$
    – Alex Trounev
    Commented Nov 28, 2021 at 12:39
  • $\begingroup$ To extend iterative solution up to tmax=3 we can put nmax=25 and use FDM solution as starting point. $\endgroup$
    – Alex Trounev
    Commented Nov 28, 2021 at 17:44
  • $\begingroup$ Also we can use as starting point X[0][t_] := Exp[-t/tmax]; Y[0][t_] := 0; $\endgroup$
    – Alex Trounev
    Commented Nov 28, 2021 at 17:54

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.