5
$\begingroup$

Bug introduced in 12.0 and fixed in 12.1

The following are the dimensionless equations of a one-step random walk with an "impure" reflecting boundary assuming only four possible states, i.e., n=0, 1, 2, 3, 4.

Equations On plugging in the equations with suitable initial conditions into Mathematica and using DSolve to solve the equations: 

  sol = DSolve[{(p0'[t] == -2*p0[t] + p1[t]), (p1'[t] == 
      2*p0[t] + p2[t] - 2*p1[t]),  (p2'[t] == 
      p1[t] + p3[t] - 2*p2[t]), p0[0] == 0, p1[0] == 0, p2[0] == 1, 
    p3[0] == 0, p0[t] + p1[t] + p2[t] + p3[t] == 1}, {p0[t], p1[t], 
    p2[t], p3[t]}, {t}];

   Plot[Evaluate[{p0[t], p1[t], p2[t], p3[t]} /. sol], {t, 0, 10}, 
 PlotLegends -> {"\!\(\*SubscriptBox[\(P\), \(0\)]\)(t)", 
   "\!\(\*SubscriptBox[\(P\), \(1\)]\)(t)", 
   "\!\(\*SubscriptBox[\(P\), \(2\)]\)(t)", 
   "\!\(\*SubscriptBox[\(P\), \(3\)]\)(t)"}, 
 PlotRange -> {{0, 10}, {0, 1}}]

It took a while but I got the plots. Then, I tried doing the same with NDSolve,

numericirb = 
  NDSolve[{(p0'[t] == -2*p0[t] + p1[t]), (p1'[t] == 
      2*p0[t] + p2[t] - 2*p1[t]),  (p2'[t] == 
      p1[t] + p3[t] - 2*p2[t]), p0[0] == 0, p1[0] == 0, p2[0] == 1, 
    p3[0] == 0, p0[t] + p1[t] + p2[t] + p3[t] == 1}, {p0[t], p1[t], 
    p2[t], p3[t]}, {t, 0, 10}];

Plot[{p0[t] /. numericirb, p1[t] /. numericirb, p2[t] /. numericirb, 
  p3[t] /. numericirb}, {t, 0, 10}, 
 PlotLegends -> {"\!\(\*SubscriptBox[\(P\), \(0\)]\)(t)", 
   "\!\(\*SubscriptBox[\(P\), \(1\)]\)(t)", 
   "\!\(\*SubscriptBox[\(P\), \(2\)]\)(t)", 
   "\!\(\*SubscriptBox[\(P\), \(3\)]\)(t)"}, 
 PlotRange -> {{0, 10}, {0, 1}}]

Here are the plots from both the outputs:

Plot from DSolve

Plot from NDSolve

I do not understand why is there a discrepancy since the equations plugged in are identical? What am I missing here? I am using version 12.0.0 for Mac OS X x86 (64-bit) on Mac OS Catalina version 10.15.4

Improve this question
$\endgroup$
4
  • $\begingroup$ In v9.0.1, the plots are the same as your second plot. Which version are you in? $\endgroup$
    – xzczd
    Apr 14, 2020 at 4:32
  • $\begingroup$ 12.0.0 for Mac OS X x86 (64-bit). Also I'm sorry about the inefficient code used for plotting. Is there any error there? $\endgroup$
    – Kk85
    Apr 15, 2020 at 16:57
  • $\begingroup$ It's a bug. You can substitute back the DSolve solution and see that the equations are not satisfied. $\endgroup$
    – Szabolcs
    Apr 15, 2020 at 17:31
  • $\begingroup$ Ah, my bad! Thank you. $\endgroup$
    – Kk85
    Apr 15, 2020 at 17:33

1 Answer 1

Reset to default
4
$\begingroup$

tl;dr This is a bug in version 12.0.0.

When there is a discrepancy between the symbolic and the numeric result, the symbolic one is usually (not always) incorrect. To make sure, the result can always be substituted back to verify it.

Here's a simpler and more concise way to write the problem:

eqn = {
   p0'[t] == -2*p0[t] + p1[t],
   p1'[t] == 2*p0[t] + p2[t] - 2*p1[t],
   p2'[t] == p1[t] + p3[t] - 2*p2[t],
   p0[0] == 0, p1[0] == 0, p2[0] == 1, p3[0] == 0,
   p0[t] + p1[t] + p2[t] + p3[t] == 1
   };

funs = {p0, p1, p2, p3};

sol = DSolve[eqn, funs, t];

nsol = NDSolve[eqn, funs, {t, 0, 10}];

Plotting:

Plot[{p0[t], p1[t], p2[t], p3[t]} /. sol // Evaluate, {t, 0, 10},
 PlotLegends -> funs, PlotRange -> All] (* symbolic *)

Plot[{p0[t], p1[t], p2[t], p3[t]} /. nsol // Evaluate, {t, 0, 10},
 PlotLegends -> funs, PlotRange -> All] (* numeric *)

Let's verify the symbolic solution by back-substitution.

Table[
  (Subtract @@@ eqn) /. First[sol] /. t -> tt,
  {tt, 0., 10}
  ] // Chop

(*

{{0, 1., -3., 0, 0, 0, 0, 0}, {0, 0.367879, -1.73576, 0, 0, 0, 0, 
  0}, {0, 0.135335, -1.27067, 0, 0, 0, 0, 0}, {0, 0.0497871, -1.09957,
   0, 0, 0, 0, 0}, {0, 0.0183156, -1.03663, 0, 0, 0, 0, 0}, {0, 
  0.00673795, -1.01348, 0, 0, 0, 0, 0}, {0, 0.00247875, -1.00496, 0, 
  0, 0, 0, 0}, {0, 0.000911882, -1.00182, 0, 0, 0, 0, 0}, {0, 
  0.000335463, -1.00067, 0, 0, 0, 0, 0}, {0, 0.00012341, -1.00025, 0, 
  0, 0, 0, 0}, {0, 0.0000453999, -1.00009, 0, 0, 0, 0, 0}}

*)

If the solution is correct, this should have yielded all zeros.

It turns out that the solution is correct and simple in version 11.3. It is also correct in version 12.1.0, but it is not simple at all. In fact, the result I get is so large that I needed to use the options MaxRecursion -> 0, PlotPoints -> 30 be able to plot it in a reasonable amount of time.

For reference, the solution returned by 11.3 is:

{{p0 -> Function[{t}, (1/
    4950967)(707281 + 
      67355776 RootSum[7 + 14 #1 + 7 #1^2 + #1^3 &, 
        E^(t #1)/(14 + 14 #1 + 3 #1^2) &] + 
      34933554 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (2 E^(t #1) + E^(t #1) #1)/(
         14 + 14 #1 + 3 #1^2) &] - 
      45308334 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (3 E^(t #1) + E^(t #1) #1)/(
         14 + 14 #1 + 3 #1^2) &] + 
      12615345 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (5 E^(t #1) + 2 E^(t #1) #1)/(
         14 + 14 #1 + 3 #1^2) &] + 
      16270472 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (
         2 E^(t #1) + 4 E^(t #1) #1 + E^(t #1) #1^2)/(
         14 + 14 #1 + 3 #1^2) &] - 
      16977753 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (
         6 E^(t #1) + 5 E^(t #1) #1 + E^(t #1) #1^2)/(
         14 + 14 #1 + 3 #1^2) &])], 
  p1 -> Function[{t}, (1/
    4950967)(1414562 + 
      83950995 RootSum[7 + 14 #1 + 7 #1^2 + #1^3 &, 
        E^(t #1)/(14 + 14 #1 + 3 #1^2) &] + 
      46111767 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (2 E^(t #1) + E^(t #1) #1)/(
         14 + 14 #1 + 3 #1^2) &] - 
      53965202 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (3 E^(t #1) + E^(t #1) #1)/(
         14 + 14 #1 + 3 #1^2) &] + 
      16881210 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (5 E^(t #1) + 2 E^(t #1) #1)/(
         14 + 14 #1 + 3 #1^2) &] + 
      23787142 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (
         2 E^(t #1) + 4 E^(t #1) #1 + E^(t #1) #1^2)/(
         14 + 14 #1 + 3 #1^2) &] - 
      25201704 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (
         6 E^(t #1) + 5 E^(t #1) #1 + E^(t #1) #1^2)/(
         14 + 14 #1 + 3 #1^2) &])], 
  p2 -> Function[{t}, (1/
    4950967)(1414562 - 
      26090495 RootSum[7 + 14 #1 + 7 #1^2 + #1^3 &, 
        E^(t #1)/(14 + 14 #1 + 3 #1^2) &] - 
      15412277 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (2 E^(t #1) + E^(t #1) #1)/(
         14 + 14 #1 + 3 #1^2) &] + 
      20106239 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (3 E^(t #1) + E^(t #1) #1)/(
         14 + 14 #1 + 3 #1^2) &] - 
      5124060 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (5 E^(t #1) + 2 E^(t #1) #1)/(
         14 + 14 #1 + 3 #1^2) &] - 
      2725034 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (
         2 E^(t #1) + 4 E^(t #1) #1 + E^(t #1) #1^2)/(
         14 + 14 #1 + 3 #1^2) &] + 
      6261439 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (
         6 E^(t #1) + 5 E^(t #1) #1 + E^(t #1) #1^2)/(
         14 + 14 #1 + 3 #1^2) &])], 
  p3 -> Function[{t}, (1/
    4950967)(1414562 - 
      125216276 RootSum[7 + 14 #1 + 7 #1^2 + #1^3 &, 
        E^(t #1)/(14 + 14 #1 + 3 #1^2) &] - 
      65633044 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (2 E^(t #1) + E^(t #1) #1)/(
         14 + 14 #1 + 3 #1^2) &] + 
      79167297 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (3 E^(t #1) + E^(t #1) #1)/(
         14 + 14 #1 + 3 #1^2) &] - 
      24372495 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (5 E^(t #1) + 2 E^(t #1) #1)/(
         14 + 14 #1 + 3 #1^2) &] - 
      37332580 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (
         2 E^(t #1) + 4 E^(t #1) #1 + E^(t #1) #1^2)/(
         14 + 14 #1 + 3 #1^2) &] + 
      35918018 RootSum[
        7 + 14 #1 + 7 #1^2 + #1^3 &, (
         6 E^(t #1) + 5 E^(t #1) #1 + E^(t #1) #1^2)/(
         14 + 14 #1 + 3 #1^2) &])]}}
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.