DSolve apparently returns wrong symbolic result

Question

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:

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?

Answer 1

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])}} *)

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])]}}*)