MathGroup Archive 2010

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

Search the Archive

Solving ODE with discontinuous shock perturbation [was: Modification

  • To: mathgroup at smc.vnet.net
  • Subject: [mg108030] Solving ODE with discontinuous shock perturbation [was: Modification
  • From: Shawn Garbett <shawn.garbett at gmail.com>
  • Date: Sun, 7 Mar 2010 04:01:01 -0500 (EST)

Here's my final solution and a new question. This solves a basic
oscillator coupled to a growth equation, with the addition of a
discontinuity, that when F reaches 5 from below, it drops to half it's
value. Using NestList NDSolve is restarted at each discontinuity, and
a method of piecewise evaluating the solution from the resulting set
of interpolating functions is given.

dA = B[t];
dB = (1 - A[t]^2) B[t] - A[t];
dF = B[t] + 0.2;

tmax = 100;

initial = {0, {A -> Function[t, 1], B -> Function[t, 1],
    F -> Function[t, 1]}};

# Note the shock is in the initial conditions for F, and the step
allows restarting from the last position.
step[x_] := With[
  {div = First[x], sol = Last[x]},
  Block[
   {interval =
     Reap[
      NDSolve[
       {A'[t] == dA,
        B'[t] == dB,
        F'[t] == dF,
        A[div] == (A[div] /. sol),
        B[div] == (B[div] /. sol),
        F[div] == (F[div] /. sol)/2},
       {A, B, F}, {t, div, tmax},
       Method -> {
         "EventLocator",
         "Event" -> F[t] - 5,
         "Direction" -> 1,
         "EventAction" :> Sow[t]
         }
       ]]},
   {interval[[2]][[1]][[1]], interval[[1]][[1]]}
   ]]

result = Rest[NestList[step, initial, 7]];

sol[t_] := Last[First[Select[result, t < First[#] &]]]

disp = Plot[#[t] /. sol[t], {t, 1, 100}, PlotRange -> All] &;
GraphicsGrid[Map[disp, {{A}, {B}, {F}}, {2}]]

##
This example works. However, when I applied it to the full problem (18
coupled equations), the EventLocator is giving me a point outside the
requested solution range.

This solves from 0 to 100 (much bigger problem, not example above)

x1 = step[initial]
{7.51929, {A -> InterpolatingFunction[{{0.,100.}},<>], ...}}

Note that the first number returned is the crossing point of the event
that triggers the discontinuity. So using this result and calling the
step function again.

x2 = step[x1]
{6.22569, {A -> InterpolatingFunction[{{7.51929, 100.0}}, <>] , ...}}

The first number (6.22569) should be inside the range (7.51929, 100.0)
of the interpolating function, but it's not. It does this consistently
for the above toy example. This was the requested solution range sent
to NDSolve at this step. So the EventLocator is returning a point
*outside* the requested solution range. I can't figure out why this is
happening. I'm guessing that there is a sensitive computation that has
a really low slope somewhere and it shoots off the range. Is there
some parameter that will constrain it to the given range? Are there
variables lurking from a previous execute of step?


  • Prev by Date: Re: learning calculus through mathematica
  • Next by Date: Re: removing non-numeric elements from the table
  • Previous by thread: Re: Dynamic function application problem
  • Next by thread: Forcing mathematica to output a certain form