Re: NDSolve - how to bypass safety chceck?
- To: mathgroup at smc.vnet.net
- Subject: [mg111394] Re: NDSolve - how to bypass safety chceck?
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Fri, 30 Jul 2010 06:57:09 -0400 (EDT)
slawek wrote:
> U=BFytkownik "Daniel Lichtblau" <danl at wolfram.com> napisa=B3 w wiadomo=B6ci grup
> dyskusyjnych:i2ok3i$7ni$1 at smc.vnet.net...
>> You neglected to mention what happened when you tried it. Was it
>> something bad (e.g. program hang or frozen computer)?
>>
>> Daniel Lichtblau
>> Wolfram Research
>
>
> I solve sets of DDE equations. The simplest example is:
>
> {x1'[t] == x1[t] - x1[t - 2 d] * x2[t - d]^2 - 2 * x1[x - d] *
> x2[x - 2 d] * x2[x - d],
> x2'[t] == x2[t] - x2[t - 2 d] * x1[t - d]^2 - 2 * x2[x - d] *
> x1[x - 2 d] * x1[x - d]}
>
> Above equations are related to integro-differential equations, which are can
> not be solved non-numerically and which are nonlinear plasma physics
> equations. The history function is a main trick to solve DDE, because it
> extended domain of time t. (This "trick" and DDE have been never applied in
> published works on the topic.)
>
> The ODE-INT "old method" take hours to obtain a solution (Fortran), the DDE
> "new method" gives a solution after about half a second (Mathematica 7, the
> same PC). Wow!
>
> There are - and it is ok - some subdomains of the time shift d, where
> solutions are chaotic. I use an "module", to generate plots and animations
> (or use Slider/Dynamic):
>
> plot[d_] := Module[{...}, NDSolve[...]; ...; Plot[...]]
> plots = Table[plot[d],{d,0.05,5.0}]
> ListAnimate[plots]
>
> For d = 5.0 Mathematica refuses give a solutions and warnig is generate d. It
> is ok, because the solution is chaotic. Nevertheless I would like have even
> completly noisy "wrong" solution instead a blank plot. Thus the question how
> to weak "quality checking" in NDSolve and to force NDSolve to always produce
> an output, even if the covergence is poor.
>
> NDSolve return when it finish computations, thus - I suppose - Infinity as
> the number of steps would generate infinite computation time and no results.
> I do not look for more accurate results. I look for an option to switch
> NDSolve to a "dumb mode".
>
> TIA
> slawek
It is pretty difficult to give a useful response in the absence of
working code. Below are modifications of your partial code that seem to
"work", but do not readily exhibit the chaos or early termination of the
solver.
You will notice that I changed "x" to "t" in the DDE independent
variables, so that it would not be a partial delayed DE; this is in
accord with your notation, where derivatives are only with respect to t.
In order to convince NDSolve to actually do anything, I also had to set
up initial conditions. And sacrifice a small reptile.
eq = {x1'[t] ==
x1[t] - x1[t - 2 d]*x2[t - d]^2 -
2*x1[t - d]*x2[t - 2 d]*x2[t - d],
x2'[t] ==
x2[t] - x2[t - 2 d]*x1[t - d]^2 -
2*x2[t - d]*x1[t - 2 d]*x1[t - d], x1[t /; t <= 0] == -2,
x2[t /; t <= 0] == 1};
plot[p_] :=
Module[{sol}, sol = NDSolve[eq /. d -> p, {x1[t], x2[t]}, {t, 0, 1}];
Plot[{x1[t], x2[t]} /. First[sol], {t, 0, 1}]]
plot[.5]
Okay, so what to do when you actually run into an early termination?
Here are a few possibilities.
(1) Use "EventHandler" as a method, and code a restart.
(2) Jettison the DDE approach, and use recurrences to approximate the
behavior. This is a way of emulating a blind fixed step approach (that
is to say, there will be no error assessment, so the noise is allowed to
stay.
(3) Possibly lowering Accuracy?PrecisionGoal might help.
(4) Possibly you can use Min and Max to restrict the sizes of
derivatives. But I'm not sure that is what matters for the situation you
describe (again, without your working example, it is quite difficult to
test anything).
Daniel Lichtblau
Wolfram Research