6
$\begingroup$

Bug introduced in 9.0 and persisting through 11.0.1 or later


I am trying to get Mathematica to calculate the correct solution $x(t)$ of the following ODE:

$\dot x = a x(t)+ct^n$

for $a>0$, $n\ge 1$, $c>0$ and $x(0) = x_0$.

My final goal is to calculate the asymptotic behavior for $t\rightarrow \infty$ (i.e. the fastest growing term). I can do so on a piece of paper and I get

$\left(x_0+\frac{c n!}{a^{n+1}}\right)e^{at}$.

I am trying to do the same thing using Mathematica with this code:

eqn = D[x[t], t] == a*x[t] + c*t^n;
sol = DSolve[{eqn, x[0] == x0}, x, t]
Simplify[eqn /. sol] (*1*)
solX[t_] = (x[t] /. sol[[1]]) (*2*)
asymptX[t_] = 
  Limit[solX[t]/(Exp[a t]), t -> \[Infinity], 
  Assumptions -> {a > 0, n >= 1, c > 0}]*Exp[a t] (*3*)

The line marked "(*1*)" returns {True}, so Mathematica thinks its solution is correct. The full solution (*2*) is given as

$\frac{e^{a t} (a t)^{-n} \left(-c t^n \Gamma (n+1,a t)+a c \Gamma (n+1,0) (a t)^n+a x_0 (a t)^n\right)}{a}$

and the asymptotic solution (*3*) that is returned is

$e^{a t} (c \Gamma (n+1)+x_0) = e^{a t}(cn! + x_0)$.

As you can see, the problem is that the factor $1/a^{n+1}$ is missing from the asymptotic solution compared to my by-hand calculation. I also calculated the numerical solution with Mathematica and then checked and compared all the solutions:

asymptXReal[t_]=(x0+c n!/a^(n+1))*Exp[a t]
Block[{a = 0.8,n=5,c=0.123,x0=0.123},
  nsol=NDSolve[{eqn,x[0]==x0},x,{t,0,15}];
  Plot[{Evaluate[x[t]/.nsol],solX[t],asymptXReal[t],asymptX[t]},
    {t,0,15},PlotLegends->Placed["Expressions",Below]]]

which yields the following graph:

Plot

You see that my manual derivation agrees asymptotically with the numerical solution, while Mathematica's full and asymptotic symbolic solutions are both off (but consistent with each other). Apparently, the asymptotic solution is wrong because the full solution is wrong to begin with. But why is that? I guess I am just making a stupid mistake here. Can anyboday help me find it?

$\endgroup$
  • $\begingroup$ @xzczd, thanks for the diligence in version tracking! :) $\endgroup$ – J. M. will be back soon Jun 28 '16 at 10:23
6
$\begingroup$

This turns out to be a bug since v9. v8.0.4 gives the correct result:

$Version
eqn = D[x[t], t] == a*x[t] + c*t^n;
bc = x[0] == x0;
sol = Simplify[DSolve[{eqn, bc}, x[t], t], {a > 0, n >= 1, c > 0}]
(* "8.0 for Microsoft Windows (64-bit) (October 24, 2011)" *)
(* {{x[t] -> a^(-1 - n) E^(a t) (a^(1 + n) x0 + c Gamma[1 + n] - c Gamma[1 + n, a t])}} *)

enter image description here

The analytic solution given in v9 and v10 isn't in consistent with the b.c.:

Simplify[x[0] == x0 /. sol, {a > 0, n >= 1, c > 0}]
(* {a == 1} *)

If you're in v10, a possible workaround is to add Assumptions to DSolve:

sol2 = DSolve[{eqn, x[0] == x0}, x, t, Assumptions-> {a > 0, n >= 1, c > 0}]
(* {{x -> Function[{t}, 
     a^(-1 - n) E^(a t) (a^(1 + n) x0 + c Gamma[1 + n] - c Gamma[1 + n, a t])]}} *)

Another workaround that also works in v9 is to find the general solution first and then solve for the constant:

sol3 = # /. Solve[#[[1, 1, -1]][0] == x0, C[1]][[1]] &@DSolve[eqn, x, t];

Simplify[sol3 /. Function -> func, {a > 0, n >= 1, c > 0}] /. func -> Function
(*{{x -> Function[{t}, E^(a t) (x0 + c Gamma[1 + n] - a^(-1 - n) c Gamma[1 + n, a t])]}}*)
$\endgroup$
  • $\begingroup$ Great idea to check the b.c.! It was kind of mysterious to me how the bug could affect DSolve and Simplify at the same time to make Mathematica think its solution was correct. Thanks for trying the different versions and the workaround. I have sent a bug report to Wolfram support. $\endgroup$ – Oguk Jun 28 '16 at 7:06
  • $\begingroup$ @Oguk What was the outcome of your bug report? $\endgroup$ – bbgodfrey Aug 13 '16 at 3:53

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.