Nonlinear complex ODE system modeling modified harmonic oscillator coupled to quantum harmonic oscillator

Question

I'm using NDSolve for a system of non-linear ODEs. Here is my code

g1 = 5;
g2 = 1;
tbar = 50; 
n = 20;
(*declare matrices*)
d0 = IdentityMatrix[n];
d1 = Table[If[i == j, i, 0], {i, 0, n - 1}, {j, 0, n - 1}];
d2 = Table[If[i == j - 1, Sqrt[j], 0], {i, 0, n - 1}, {j, 0, n - 1}];
d3 = Table[If[i == j + 1, Sqrt[j + 1], 0], {i, 0, n - 1}, {j, 0, n - 1}];

(*declare initial variables*)
{x0, p0} = {0, 0};
initc = a[0] == (x0 + I*p0)/Sqrt[2];
{expx0, expp0} = {1, 1};
\[Alpha] = (expx0 + I*expp0)/Sqrt[2];
initq = c[0] == Chop[Normalize@Table[\[Alpha]^i/Sqrt[i!], {i, 0, n - 1}]];

(*differential equations*)

sol = NDSolve[{a'[t] == -I (g1*a[t] + g2/2 c[t].(d2 + d3).Conjugate[c[t]]), 
    c'[t] == -I*c[t].(d1 + 0.5*d0 + g2*(Re@a[t]) (d2 + d3)), 
    initc, initq}, {a, c}, {t, 0, tbar}];

(*plotting*)
asol[t_] := a[t] /. sol[[1]]
csol[t_] := c[t] /. sol[[1]]
Plot[{Re@asol[t], Im@asol[t]}, {t, 0, tbar}, PlotRange -> All]

I want to be able to get reliable results for this system for a range of values of the parameters $g1$ and $g2$. The problem arises when I use values which make the non-linear term (the second term) in the first equation large. For example, when I use $g1=5, g2=1$, the solutions are good but when I use $g1=0.5, g2=1$, I can already see that there is something wrong (the values are too large). So essentially the problem seems to be in the ratio $g2/g1$ and I suspect this is due to the non-linearity and perhaps something goes wrong with the NDSolve algorithm (the run times are also much longer for the problematic cases). Is there a way I can make this code work for a larger parameter range?

EDIT: what do I mean by "too large"? I mean the values of Re@asol and Im@asol being large, to see why this is an issue I need some context for this system. This is a problem from physics of a modified harmonic oscillator (MHO) coupled to a quantum harmonic oscillator (QHO). Energy conservation should hold for this system i.e. MHO+QHO+interaction is a constant.This can be checked with the following code

e1[t_] := g1*(Abs@asol[t])^2;
e2[t_] := Chop[csol[t].(d1 + 0.5*d0).Conjugate[csol[t]]];
eint[t_] := 
  Chop[csol[t].(g2*(Re@asol[t])*(d2 + d3)).Conjugate[
    csol[t]]];
etotal[t_] := e1[t] + e2[t] + eint[t];
{etotal[0], etotal[tbar]}

Energy conservation begins to fail as I go to g2=0.5 (try g2=0.1 for a more dramatic demonstration of this).

Answer 1

First of all, I'd like to emphasize again, the initial value problem (IVP) solver for ODE built in NDSolve is very robust and should always be the last thing to adjust, so please make sure the equation system itself is correct.

Anyway, let alone the correctness of system, if you care about the invariant, consider Projection method (The following takes about 350 seconds on my laptop for g1 = 0.1; g2 = 1;):

invar = g1*(Abs@a[t])^2 + c[t] . (d1 + 0.5*d0) . Conjugate[c[t]] + 
  c[t] . (g2*(Re@a[t])*(d2 + d3)) . Conjugate[c[t]]
solinvar = 
 NDSolve[{a'[t] == -I (g1*a[t] + g2/2 c[t] . (d2 + d3) . Conjugate[c[t]]), 
   c'[t] == -I*c[t] . (d1 + 0.5*d0 + g2*(Re@a[t]) (d2 + d3)), initc, initq}, {a, 
   c}, {t, 0, tbar}, 
  Method -> {"Projection", Method -> "ImplicitRungeKutta", "Invariants" -> invar}, 
  MaxSteps -> Infinity]

Plot[invar /. sol // Evaluate, {t, 0, tbar}]

Still, it's not perfectly conserved, but much better than the default. For comparison:

solref = NDSolve[{a'[t] == -I (g1*a[t] + g2/2 c[t] . (d2 + d3) . Conjugate[c[t]]), 
     c'[t] == -I*c[t] . (d1 + 0.5*d0 + g2*(Re@a[t]) (d2 + d3)), initc, initq}, {a, 
     c}, {t, 0, tbar}, MaxSteps -> Infinity]; // AbsoluteTiming

(* {33.922384, Null} *)

Plot[invar /. solref // Evaluate, {t, 0, tbar}]

However, I doubt if this really makes sense, given that this solution is virtually the same as the one given by default method:

Plot[ReIm@a[t] /. {solinvar[[1]], solref[[1]]} // Evaluate, {t, 0, tbar}, 
 PlotRange -> All, PlotStyle -> {Automatic, Automatic, Dashed, Dashed}]