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

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