# Solving simultaneous differential equations using eigen value method

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]

• I don't see what the issue is to be honest. Have you tried the obvious? Eigenvalues[Smatrix] and Eigenvectors[Smatrix] with Smatrix being the matrix you wrote down? – DiSp0sablE_H3r0 May 21 at 7:52
• I have tried and right now I am getting results for N=2. But not for large N. Of course they are the values of S matrix – Jasmine May 21 at 8:08
• So, to rephrase: is the real question how to write a command to generate the S-matrix for general values of N? – DiSp0sablE_H3r0 May 21 at 8:10
• Nop. I want to solve simultaneous ODE using the method I described above – Jasmine May 21 at 8:11
• Jasmine, can you, please, include the code you have tried so far, or what exists in other related questions that you are having trouble with? It is hard to help debug and solve such problems otherwise. – CA Trevillian May 21 at 8:17

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.

• This is awesome and the behavior is expected. But in the question, I have mentioned NN=5 and tmax=3000 It's taking very long time! Any way to reduce time – Jasmine May 22 at 10:58
• It is no surprise that it takes very long if you calculate nearly a million of periods of periods. But is this really necessary? What do you want to do with it? Some background is needed. – Daniel Huber May 22 at 12:34
• I found a paper in the arXiv and trying to replicate the Fig1 in the paper, <arxiv.org/abs/1906.04900> They are using a different time scale! – Jasmine May 22 at 13:43
• Is there a possibility to rescale x axis in such a way to get a plot with new time, $$t=(\pi/4)*t’/Tn$$ : Here t’ is the old time from 0 to 1000, Tn is a particular time in which the value of B_0(t’)*conjugate[B_0(t’)] and B_N(t’)*conjugate[B_N(t’)] nearly equal to .5 – Jasmine May 22 at 14:43
• Do you want only to draw the known solution y[t'] with a new time scale t ? For this you would plot: f[4 Tn t / Pi] If you want equations in t you have to replace multiply t' and dt' by a factor of: 4 Tn / Pi – Daniel Huber May 22 at 15:44