Solving simultaneous differential equations using eigen value method

Question

I wish to solve the following set of ODE. $$i\frac{d}{dt}B_{n}\left(t\right) =f\sqrt{\left(P-n\right)\left(n+1\right)}B_{n+1}\left(t\right)+f\sqrt{n\left(P-n+1\right)}B_{n-1}\left(t\right) + Y\left[\left(P-n\right)\left(P-n-1\right)+n\left(n-1\right)\right]B_{n}\left(t\right)$$ The question is similar to the Solving n simultaneous differential equation.It's working perfectly for N=2 case. But the method is not giving any results for larger N with larger time say t=0 to 3000, and varying constants.

Values of constants are $f=20$ , $Y=334$. Also $n=0,1,2,...P$

When we consider a simple case, $P=2$ we get 3 simulatneous ODE as follows. $$\frac{d}{dt}B_{0}\left(t\right) =-iYB_{0}\left(t\right)-i\sqrt{2}fB_{1}\left(t\right)$$

$$\frac{d}{dt}B_{1}\left(t\right) =-i\sqrt{2}fB_{0}\left(t\right)-i\sqrt{2}fB_{2}\left(t\right)$$

$$\frac{d}{dt}B_{2}\left(t\right) =-i\sqrt{2}fB_{1}\left(t\right)-iYB_{2}\left(t\right)$$

The idea is to solve the above equation by finding eigen values and eigen vectors. We can convert the above equation in a matrix form as follows: $$\begin{pmatrix}\frac{d}{dt}B_{0}\left(t\right)\\ \frac{d}{dt}B_{1}\left(t\right)\\ \frac{d}{dt}B_{2}\left(t\right) \end{pmatrix}=\left(\begin{array}{ccc} -iY & -i\sqrt{2}f & 0\\ -i\sqrt{2}f & 0 & -i\sqrt{2}f\\ 0 & -i\sqrt{2}f & -iY \end{array}\right)\begin{pmatrix}B_{0}\left(t\right)\\ B_{1}\left(t\right)\\ B_{2}\left(t\right) \end{pmatrix}$$

Let $$S=\left(\begin{array}{ccc} -iY & -i\sqrt{2}f & 0\\ -i\sqrt{2}f & 0 & -i\sqrt{2}f\\ 0 & -i\sqrt{2}f & -iY \end{array}\right)$$

For solving the set of differential equations we need to find the eigen values and eigenvectors of $S$. When we calcualte them, it turns out to be 3 complex and distinct eigen values and 3 eigenvectors corresponding to each eigenvalues. Thus the solution will be:

$$B_{n}\left(t\right)=\sum_{n}G_{n}\overline{j}e^{\lambda t}$$

provided: $\overline{j}$ are eigen vectors and $\lambda$ is the corresponding eigenvalue.

$G_{n}$ can be solved using the initial condition, $$B_{n}\left(0\right)=\begin{pmatrix}1\\ 0\\ 0 \end{pmatrix}$$

The above is just $P=2$ case. Any way to extend the above for $n$ case?

After that I need to plot the modulus square of the coeficcients with respect to time in a single plot.

What I have tried using NDSolve is:

 NN = 5;

f = 20



Y = 334



 N1 = Sqrt[(NN - n)*(n + 1)];

N2 = Sqrt[n*(NN - n + 1)];
 N3 = ((NN - n)*(NN - n - 1)) + (n*(n - 1));

FT = f*N1



ST = f*N2



TT = (Y/2)*N3


 odes = Table[
   I ToExpression["M" <> ToString[n]]'[t] == 
    FT*ToExpression["M" <> ToString[n + 1]][t] + 
     ST*ToExpression["M" <> ToString[n - 1]][t] + 
     TT*ToExpression["M" <> ToString[n]][t], {n, 0, NN}];

deps = Table[ToExpression["M" <> ToString[n]][t], {n, 0, NN}]

 {M0[t], M1[t], M2[t], M3[t], M4[t], M5[t]}

 ic = {M0[0] == 1, M1[0] == 0, M2[0] == 0, M3[0] == 0, M4[0] == 0, M5[0] == 0}


MH = NDSolve[{odes, ic}, deps, {t, 0, 3000}] 
  


Plot[{Evaluate[(M0[t]*Conjugate[M0[t]]) /. MH], 
  Evaluate[(M1[t]*Conjugate[M1[t]]) /. MH], 
  Evaluate[(M2[t]*Conjugate[M2[t]]) /. MH], 
  Evaluate[(M3[t]*Conjugate[M3[t]]) /. MH], 
  Evaluate[(M4[t]*Conjugate[M4[t]]) /. MH], 
  Evaluate[(M5[t]*Conjugate[M5[t]]) /. MH]}, {t, 0, 3000}, PlotRange -> All]

Answer 1

I think your problem is that you have very fast oscillating solutions besides very slow varying ones.

I calculated (with some simplifications) in reasonable time the solution for NN=4 and tmax=1000:

NN = 4;
tmax = 1000;
f = 20;
Y = 334;
N1 = Sqrt[(NN - n)*(n + 1)];
N2 = Sqrt[n*(NN - n + 1)];
N3 = ((NN - n)*(NN - n - 1)) + (n*(n - 1));
FT = f*N1;
ST = f*N2;
TT = (Y/2)*N3;

odes = Table[
   I ToExpression["M" <> ToString[n]]'[t] == 
    FT*ToExpression["M" <> ToString[n + 1]][t] + 
     ST*ToExpression["M" <> ToString[n - 1]][t] + 
     TT*ToExpression["M" <> ToString[n]][t], {n, 0, NN}];

deps = Table[ToExpression["M" <> ToString[n]][t], {n, 0, NN}];
ic = Prepend[Thread[Rest[deps /. t -> 0] == 0], M0[0] == 1];

(MH = NDSolve[Flatten@{odes, ic}, deps, {t, 0, tmax}][[1]] ;
  
  Plot[#, {t, 0, tmax}, 
     PlotRange -> All] & /@ ((deps Conjugate[deps]) /. MH)
  ) // Timing

As you can see, it took 138 sec and the main point is, that solutions 2,3,4 oscillate very fast. How really fast you can see, if we plot solution 2 for 1/1000 of tmax, time t=0 to t=1:

Plot[Evaluate[(deps Conjugate[deps]) /. MH[[2]]], {t, 0, 1}, PlotRange -> All]

This should illuminate the problem enough.

What can be done? Do you need the accurate solutions of the fast varying parts? If yes, you could calculate the fast changing parts over a small time interval and replace the slow moving parts by constants. On the other hand, the slow moving parts seem to feel only an average of the fast parts, so you could replace the fast parts in the equations for the slow parts by an average.