MathGroup Archive 2008

[Date Index] [Thread Index] [Author Index]

Search the Archive

How to do efficient NDSolve Error control

  • To: mathgroup at smc.vnet.net
  • Subject: [mg91757] How to do efficient NDSolve Error control
  • From: pratip <pratip.official at gmail.com>
  • Date: Sun, 7 Sep 2008 05:37:03 -0400 (EDT)

Lets take this system of ODE and solve with DSolve

exsol = DSolve[{x1'[t] == x2[t],
x2'[t] == c^2 x1[t] - (c^2 + p^2) Sin[p t], x1[0] == 0,
x2[0] == Pi}, {x1[t], x2[t]}, t] // FullSimplify

solution is as follows
{{x1[t] -> Sin[p t] + ((-p + \[Pi]) Sinh[c t])/c,
x2[t] -> p Cos[p t] + (-p + \[Pi]) Cosh[c t]}}

Here is the values of the parameters

dp = {c -> 100, p -> 2}

Now solve with NDSolve

sol = NDSolve[{x1'[t] == x2[t],
x2'[t] == c^2 x1[t] - (c^2 + p^2) Sin[p t], x1[0] == 0,
x2[0] == Pi} /. dp, {x1, x2}, {t, 0, 2}]

Lets Plot the residual error
Res = {-x1'[t] + x2[t], -x2'[t] + c^2 x1[t] - (c^2 + p^2) Sin[p t]};
Plot[Evaluate[Res /. sol /. dp], {t, 0, 2}]

One can see the exponential increase in error. Please tell me if it is possible to manage this error for any parameter values. I want to know how I can use NDSolve to simulate such ODE with very little residual error just like any other normal system of ODE. Should I choose some special method options?

Another very funny thing

der = D[exsol, t];
Res /. exsol /. der /. dp // Flatten // FullSimplify

gives {0,0}. But lets take the residual with out simplifying.
res = Res /. exsol /. der /. dp // Flatten;
Plot[Evaluate[res], {t, 0, 1.2}, PlotStyle -> Thick]

The Plot shows how badly exponential error shows up in the boundary. Is it a problem with Plot command? why it can not understand without simplifying that "res={0,0}". I hope you guys have an answer.


  • Prev by Date: Re: LU Decomposition w/o Pivoting
  • Next by Date: Re: phase-space versus controlling parameter surface
  • Previous by thread: Re: Use of delayed assignment, :=, with a list
  • Next by thread: Re: Mathematica and F#