1
$\begingroup$

I am solving 6 differential equations in matrix form as:

\[Theta] = 34.5*\[Pi]/180;
m2 = 1/2*7.3*10^-5*{{Cos[
     2 \[Theta]], -Sin[2 \[Theta]]}, {-Sin[2 \[Theta]], -Cos[
      2 \[Theta]]}}
M = Array[Subscript[Subscript[\[Rho], 1], #1, #2][t] &, {2, 2}]
M2 = Array[Subscript[c, #1, #2][t] &, {2, 2}]
M0 = {{1, 0.00001}, {0.00001, 0}};
ci = Thread[Flatten[M] == Flatten[M0]] /. {t -> 0};
M20 = {{2, 0.00001}, {0.00001, 0}};
di = Thread[Flatten[M2] == Flatten[M20]] /. {t -> 0};

s = NDSolve[{ 


 I D[M, t] == (m2/10 + 10^5*2*M2).M - M.(m2/10 + 10^5*2*M2), 
   I D[M2, t] == (m2/11 + 10^5*2*M).M2 - M2.(m2/11 + 10^5*2*M)}, ci, 
  di, {Variables[M], Variables[M2]}, {t, 0, 10^7}]

But it is showing the error as :

Subscript[c, 1,1][0]==2 cannot be used as a variable. >> Can someone please help me with this?

$\endgroup$
  • 2
    $\begingroup$ Your initial conditions should be enclosed within your curly braces with the equations. I reached the maximum number of steps at t=0.096, so you will have to investigate that. $\endgroup$ – Tim Laska Jan 9 at 11:40
  • $\begingroup$ Thank you. It worked. $\endgroup$ – Madhurima Chakraborty Jan 10 at 5:45
4
$\begingroup$

We can normalize the equations by 10^5. So that the oscillation period is of the order of 1, we make the change t1->t*10^5. Then it is possible to integrate up to 10^5, since the dynamics are also visible on this scale.

\[Theta] = 34.5*\[Pi]/180;
m2 = 1/2*7.3*10^-5*{{Cos[
      2 \[Theta]], -Sin[2 \[Theta]]}, {-Sin[2 \[Theta]], -Cos[
       2 \[Theta]]}};
M = Array[Subscript[Subscript[\[Rho], 1], #1, #2][t] &, {2, 2}];
M2 = Array[Subscript[c, #1, #2][t] &, {2, 2}];
M0 = {{1, 0.00001}, {0.00001, 0}};
ci = Thread[Flatten[M] == Flatten[M0]] /. {t -> 0};
M20 = {{2, 0.00001}, {0.00001, 0}};
di = Thread[Flatten[M2] == Flatten[M20]] /. {t -> 0};

s = NDSolve[{I D[M, t] == (m2/10^6 + 2*M2).M - M.(m2/10^6 + 2*M2), 
   I D[M2, t] == (m2/11/10^5 + 2*M).M2 - M2.(m2/11/10^5 + 2*M), ci, 
   di}, {Variables[M], Variables[M2]}, {t, 0, 10^5}, 
  MaxSteps -> Infinity]

Two-scale visualization

    {Plot[Evaluate[Variables[M] /. First[s][[1]] // Im], {t, 0, 10^1}, 
  PlotRange -> All, Frame -> True, Axes -> False, 
  PlotLegends -> Variables[M]], 
 Plot[Evaluate[Variables[M2] /. First[s][[2]] // Im], {t, 0, 10^1}, 
  PlotRange -> All, Frame -> True, Axes -> False, 
  PlotLegends -> Variables[M2]]}

{Plot[Variables[M] /. First[s][[1]] // Im, {t, 0, 10^5}, 
  PlotRange -> All, Frame -> True, Axes -> False], 
 Plot[Variables[M2] /. First[s][[2]] // Im, {t, 0, 10^5}, 
  PlotRange -> All, Frame -> True, Axes -> False]}

Figure 1

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you. What I need to do to plot the (2,2) element of the matrix M? $\endgroup$ – Madhurima Chakraborty Jan 10 at 5:44
  • $\begingroup$ Use s = NDSolve[{I D[M, t] == (m2/10^6 + 2*M2).M - M.(m2/10^6 + 2*M2),I D[M2, t] == (m2/11/10^5 + 2*M).M2 - M2.(m2/11/10^5 + 2*M), ci, di}, Flatten[{Variables[M], Variables[M2]}], {t, 0, 10^5}, MaxSteps -> Infinity] and then {Plot[Variables[M][[4]] /. First[s] // Re, {t, 0, 10^5}], Plot[Variables[M][[4]] /. First[s] // Im, {t, 0, 10^5}]} $\endgroup$ – Alex Trounev Jan 10 at 11:24

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.