在解多耦合常微分方程过程中，刚性系统出现奇点

Question

我在解一个这样的耦合系统，因为系统参数数量级差距比较大，因此系统是刚性的，
Clear["Global`*"]
ga = 559.*10^16;
m = 5.*10^-11;
\[Gamma]m = 24.*10^-8;
gc = 559.*10^16;
\[Omega]m = 735.*10^5;
\[Omega]c = 193.*10^7;
\[Omega]a = 100.*10^12;
\[CapitalDelta]a = \[Omega]m;
\[CapitalDelta]c = 0.;
\[Phi] = 0.;
\[Kappa]a = 2.*\[Omega]m/5.;
\[Kappa]c = 4.*\[Omega]m/5.;
hbar = 663.*10^-36/(2*\[Pi]);
sol = NDSolve[{q'[t] == p[t]/m,
   p'[t] == -\[Gamma]m p[t] + hbar ga (ar[t]^2 + ai[t]^2) + 
     hbar gc (cr[t]^2 + ci[t]^2),
   ar'[t] == -\[Kappa]a ar[t] + (\[CapitalDelta]a - ga q[t]) ai[
       t] + \[Omega]a,
   ai'[t] == -\[Kappa]a ai[t] - (\[CapitalDelta]a - ga q[t]) ar[t],
   cr'[t] == -\[Kappa]c cr[t] + (\[CapitalDelta]c - gc q[t]) ci[
       t] + \[Omega]c Cos[\[Phi]],
   ci'[t] == -\[Kappa]c ci[t] - (\[CapitalDelta]c - gc q[t]) cr[
       t] - \[Omega]c Sin[\[Phi]],
   \[Epsilon]q'[t] == \[Epsilon]p[t]/m,
   \[Epsilon]p'[t] == -\[Gamma]m \[Epsilon]p[t] - m \[Omega]m^2 + 
     2 hbar ga (ar[t] \[Epsilon]ar[t] + ai[t] \[Epsilon]ai[t]) + 
     2 hbar gc (cr[t] \[Epsilon]cr[t] + ci[t] \[Epsilon]ci[t]),
   \[Epsilon]ar'[t] == -\[Kappa]a \[Epsilon]ar[t] - 
     ga ai[t] \[Epsilon]q[
       t] + (\[CapitalDelta]a - ga q[t]) \[Epsilon]ai[t],
   \[Epsilon]ai'[t] == -\[Kappa]a \[Epsilon]ai[t] + 
     ga ar[t] \[Epsilon]q[
       t] - (\[CapitalDelta]a - ga q[t]) \[Epsilon]ar[t],
   \[Epsilon]cr'[t] == -\[Kappa]c \[Epsilon]cr[t] - 
     gc ci[t] \[Epsilon]q[
       t] + (\[CapitalDelta]c - gc q[t]) \[Epsilon]ci[t],
   \[Epsilon]ci'[t] == -\[Kappa]c \[Epsilon]ci[t] + 
     gc cr[t] \[Epsilon]q[
       t] - (\[CapitalDelta]c - gc q[t]) \[Epsilon]cr[t], p[0] == 0, 
   q[0] == 0, ar[0] == 0, ai[0] == 0, cr[0] == 0, 
   ci[0] == 0, \[Epsilon]p[0] == 0, \[Epsilon]q[0] == 
    0, \[Epsilon]ai[0] == 0.1, \[Epsilon]ar[0] == 
    0.1, \[Epsilon]ci[0] == 0, \[Epsilon]cr[0] == 0}, {q, p, ar, ai, 
   cr, ci, \[Epsilon]p, \[Epsilon]q, \[Epsilon]ai, \[Epsilon]ar, \
\[Epsilon]cr, \[Epsilon]ci}, {t, 0, 1*10^-6}, MaxSteps -> Infinity]
Plot[Log10[((ai[t] + \[Epsilon]ai[t])^2 + (ar[t] + \[Epsilon]ar[
        t])^2 - (ai[t]^2 + ar[t]^2))] /. sol, {t, 0, 5*10^-7}]

在解方程时出现了如下情况 

 At t == 6.502718555087617`*^-7, step size is effectively zero; singularity or stiff system suspected.

在这种情况下，我该如何处理呢？

通过查询帮助及SE上的例子 我发现用

Method -> "StiffnessSwitching"

的方法可以解决刚性矩阵和奇点这个问题

Clear["Global`*"]
Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
ga = 559.*10^16;
m = 5.*10^-11;
\[Gamma]m = 24.*10^-8;
gc = 559.*10^16;
\[Omega]m = 735.*10^5;
\[Omega]c = 193.*10^7;
\[Omega]a = 100.*10^12;
\[CapitalDelta]a = \[Omega]m;
\[CapitalDelta]c = 0.;
\[Phi] = 0.;
\[Kappa]a = 2.*\[Omega]m/5.;
\[Kappa]c = 4.*\[Omega]m/5.;
hbar = 663.*10^-36/(2*\[Pi]);
sol = NDSolve[{q'[t] == p[t]/m, 
   p'[t] == -\[Gamma]m p[t] + hbar ga (ar[t]^2 + ai[t]^2) + 
     hbar gc (cr[t]^2 + ci[t]^2), 
   ar'[t] == -\[Kappa]a ar[t] + (\[CapitalDelta]a - ga q[t]) ai[
       t] + \[Omega]a, 
   ai'[t] == -\[Kappa]a ai[t] - (\[CapitalDelta]a - ga q[t]) ar[t], 
   cr'[t] == -\[Kappa]c cr[t] + (\[CapitalDelta]c - gc q[t]) ci[
       t] + \[Omega]c Cos[\[Phi]], 
   ci'[t] == -\[Kappa]c ci[t] - (\[CapitalDelta]c - gc q[t]) cr[
       t] - \[Omega]c Sin[\[Phi]], \[Epsilon]q'[t] == \[Epsilon]p[t]/
     m, \[Epsilon]p'[t] == -\[Gamma]m \[Epsilon]p[t] - 
     m \[Omega]m^2 + 
     2 hbar ga (ar[t] \[Epsilon]ar[t] + ai[t] \[Epsilon]ai[t]) + 
     2 hbar gc (cr[t] \[Epsilon]cr[t] + 
        ci[t] \[Epsilon]ci[t]), \[Epsilon]ar'[
     t] == -\[Kappa]a \[Epsilon]ar[t] - 
     ga ai[t] \[Epsilon]q[
       t] + (\[CapitalDelta]a - ga q[t]) \[Epsilon]ai[
       t], \[Epsilon]ai'[t] == -\[Kappa]a \[Epsilon]ai[t] + 
     ga ar[t] \[Epsilon]q[
       t] - (\[CapitalDelta]a - ga q[t]) \[Epsilon]ar[
       t], \[Epsilon]cr'[t] == -\[Kappa]c \[Epsilon]cr[t] - 
     gc ci[t] \[Epsilon]q[
       t] + (\[CapitalDelta]c - gc q[t]) \[Epsilon]ci[
       t], \[Epsilon]ci'[t] == -\[Kappa]c \[Epsilon]ci[t] + 
     gc cr[t] \[Epsilon]q[
       t] - (\[CapitalDelta]c - gc q[t]) \[Epsilon]cr[t], p[0] == 0, 
   q[0] == 0, ar[0] == 0, ai[0] == 0, cr[0] == 0, 
   ci[0] == 0, \[Epsilon]p[0] == 0, \[Epsilon]q[0] == 
    0, \[Epsilon]ai[0] == 10, \[Epsilon]ar[0] == 
    10, \[Epsilon]ci[0] == 0, \[Epsilon]cr[0] == 0}, {q, p, ar, ai, 
   cr, ci, \[Epsilon]p, \[Epsilon]q, \[Epsilon]ai, \[Epsilon]ar, \
\[Epsilon]cr, \[Epsilon]ci}, {t, 0, 6*10^-6}, MaxSteps -> Infinity, 
  Method -> "StiffnessSwitching"]
Plot[Log10[((ai[t] + \[Epsilon]ai[t])^2 + (ar[t] + \[Epsilon]ar[
         t])^2 - (ai[t]^2 + ar[t]^2))] /. sol, {t, 0, 5*10^-6}]

但是在计算的后端 系统该有的震荡却消失了 

其中 我取后面极小时间内的图为 

不知道这里有没有更好的方法应付这里的这种情况

Comments

今天换了台电脑，再算就可以了，可是我把计算区间换成5*10^-6貌似又不行了。。。 这个警告会影响计算结果吗？

---

"StiffnessSwitching"不能随便用，因为当方程的解本身就含奇点时（也就是说“方程算到某个位置就算不下去了”并不是数值计算方法不当时）它依旧可能以某种方式绕过奇点继续计算最终给出错误的结果。对于这个问题我个人印象比较深的一个例子是（错误参数造成的）黑洞：万有引力方程参数比较多，因此很容易给错参数（确切地说是搞错参数的单位），这种错误可能导致引力过大，使得中心恒星实质上变成了黑洞，进而使方程的解中出现奇点，如果在此时没能及时意识到奇点出现的正确原因，而一味地在求解方法上找理由以至找上了StiffnessSwitching，你就会观察到行星在撞上黑洞时被黑洞给弹了回来——那个弹回的图像和你新贴的这个图颇有几分相似。

---

原来是这样，我看SE上用用shooting method的。我这个问题是否合适？

---

……你这又不是边值问题，不关打靶法的事儿啊。

---

是我病急乱投医了，谢谢

Answer 1

请问这个问题最后是怎么解决的？我也碰到这个问题，在自变量取值小的时候可以解不报错，但是取值大的时候就有这样的报错。就是下面的链接：

https://mmaqa.com/index.php/1831/ndsolvevalue%E5%87%BA%E7%8E%B0%E8%B6%85%E5%87%BA-9108306745948461-effectively-singularity      这里的代码如果M5[1/174]是10^33量级以下都没问题，但是一到34次方及34以上都报错，可以得到数值结果，但是都是100次方量级以上，明显不对。

 

Comments

首先这不是一个答案，所以-1，然后，我在上面的评论几乎全部可以用在你这个问题里。

---

上面的我都试过了，都不能解决，所以我觉得是不是我理解不够所以才追问的

---

“……把系数变成准确数再把WorkingPrecision调高一点试试，还不行通常就表明方程的刚性和极端数据无关，如果模型本不该有这种行为那就该考虑方程是不是列错了——对于复杂方程组这种情况非常多见。”