5
$\begingroup$

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] == 0,
    y[0] == 20
}, {x, y}, t]], {t, 0, 50}, PlotRange -> Full]

I get

DSolve

where DSolve gives the solution

{{x->Function[{t},-(1/70) E^(-50 t-7 Sqrt[51] t) (707-99 Sqrt[51]+707 E^(14 Sqrt[51] t)+99 Sqrt[51] E^(14 Sqrt[51] t)-1414 E^(14 Sqrt[51] t+(50-7 Sqrt[51]) t))],y->Function[{t},10/7 E^(-50 t-7 Sqrt[51] t) (7-Sqrt[51]+7 E^(14 Sqrt[51] t)+Sqrt[51] E^(14 Sqrt[51] 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] == 0,
    y[0] == 20
}, {x, y}, {t, 0, 50}]], {t, 0, 50}, PlotRange -> Full]

It seems like DSolve is incorrect, but why?

$\endgroup$
15
  • 1
    $\begingroup$ What if you change the 0.01 to 1/100, and set Method -> "StiffnessSwitching" in NDSolve[]? $\endgroup$ Jun 12, 2015 at 0:29
  • $\begingroup$ @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! $\endgroup$
    – user541686
    Jun 12, 2015 at 0:33
  • $\begingroup$ @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? $\endgroup$
    – user541686
    Jun 12, 2015 at 0:41
  • $\begingroup$ Confirmed. On Mma v9 both give the same result (your second plot) $\endgroup$ Jun 12, 2015 at 0:44
  • $\begingroup$ @belisarius: It's correct even when you extend it to {t,0,20}? $\endgroup$
    – user541686
    Jun 12, 2015 at 0:46

1 Answer 1

8
$\begingroup$

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] == 0, y[0] == 20}, {x, 
    y}, t];

Now consider y[9.] vs. y[9] and y[9.`20]:

y[9.] /. sol

N::meprec: Internal precision limit $MaxExtraPrecision = 50.` reached while evaluating -Sqrt[51]. >>

Mathematica graphics

y[9] /. 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]

Mathematica graphics

$\endgroup$

Your Answer

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

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