# Why does DSolve give a different result than NDSolve?

When I try DSolve to solve a system:

Plot[Evaluate[y[t] /. DSolve[{
x'[t] == y[t] / 100,
y'[t] == -100 x[t] - 100 y[t] + 2020,
x == 0,
y == 20
}, {x, y}, t]], {t, 0, 50}, PlotRange -> Full]


I get where DSolve gives the solution

{{x->Function[{t},-(1/70) E^(-50 t-7 Sqrt t) (707-99 Sqrt+707 E^(14 Sqrt t)+99 Sqrt E^(14 Sqrt t)-1414 E^(14 Sqrt t+(50-7 Sqrt) t))],y->Function[{t},10/7 E^(-50 t-7 Sqrt t) (7-Sqrt+7 E^(14 Sqrt t)+Sqrt E^(14 Sqrt t))]}}


$$x(t) = \frac{e^{-7 \sqrt{51} t-50 t}}{70} \times \left(707 e^{14 \sqrt{51} t}+99 \sqrt{51} e^{14 \sqrt{51} t}-1414 e^{\left(50-7 \sqrt{51}\right) t+14 \sqrt{51} t}+707-99 \sqrt{51}\right)$$ $$y(t) = \frac{10}{7} e^{-7 \sqrt{51} t-50 t} \left(7 e^{14 \sqrt{51} t}+\sqrt{51} e^{14 \sqrt{51} t}+7-\sqrt{51}\right)$$

which is a different plot than when I use NDSolve to solve the same system:

Plot[Evaluate[y[t] /. NDSolve[{
x'[t] == y[t] / 100,
y'[t] == -100 x[t] - 100 y[t] + 2020,
x == 0,
y == 20
}, {x, y}, {t, 0, 50}]], {t, 0, 50}, PlotRange -> Full] It seems like DSolve is incorrect, but why?

• What if you change the 0.01 to 1/100, and set Method -> "StiffnessSwitching" in NDSolve[]? – J. M. will be back soon Jun 12 '15 at 0:29
• @J. M.: Interesting, simply changing 0.01 y[t] to y[t]/100 makes DSolve match NDSolve. Seems like it turns from numeric into symbolic, so I guess that means NDSolve is actually the correct one! Thanks! – Mehrdad Jun 12 '15 at 0:33
• @J. M.: Actually, I see something even more bizarre now. When I extend the time of DSolve to {t,0,20}, I see this -- why is that? – Mehrdad Jun 12 '15 at 0:41
• Confirmed. On Mma v9 both give the same result (your second plot) – Dr. belisarius Jun 12 '15 at 0:44
• @belisarius: It's correct even when you extend it to {t,0,20}? – Mehrdad Jun 12 '15 at 0:46

It looks like classic catastrophic round-off error. (Look at those exponents on $e$!).

{sol} = DSolve[{x'[t] == y[t]/100,
y'[t] == -100 x[t] - 100 y[t] + 2020, x == 0, y == 20}, {x,
y}, t];


Now consider y[9.] vs. y and y[9.20]:

y[9.] /. sol


N::meprec: Internal precision limit \$MaxExtraPrecision = 50. reached while evaluating -Sqrt. >> y /. sol // N
(*  18.46310871964039  *)

y[9.20] /. sol
(*  18.463108719640715  *)


Try using a higher working precision:

Plot[Evaluate[y[t] /. sol], {t, 0, 50}, PlotRange -> Full,
WorkingPrecision -> 20]
` 