Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm trying to tweak the NDSolve function to perform one elementary integration step (using some explicitly selected stepping algorithm via Method option). Such a possibility is crucial for me, since after each performed step I'd like to check for some condition and modify state vector, then continue the integration. To the best of my knowledge such an operation is not possible with the 'pure' NDSolve, therefore I'm looking for an alternative solution.

So far, I managed to come up with the following solution (please note that the order of a method is selected to mimic the case with a lot of internal steps, this is equivalent to the higher order and longer integration period - tmax):

eq = {{x''[t] == -x[t], x[0] == 0, x'[0] == 1}, x, t};
tmax = 1;
method = {"ExplicitRungeKutta", "DifferenceOrder" -> 2};

Block[{state, options, timer},
  options = {MaxSteps -> Infinity, Method -> method};
  state = First@NDSolve`ProcessEquations[Sequence @@ eq, options];
  timer = First@Timing[NDSolve`Iterate[state, tmax]];
  {timer, state@"CurrentTime"["Forward"], state@"TimeStepsUsed"["Forward"], state@"SolutionVector"["Forward"]}
]

(* ==> {0.048003, 1., 7032, {0.841471, 0.540302}} *)


Block[{state, options, ctime, ic, dt, count, timer},
  options = {MaxSteps -> 1, Method -> method};
  state = First@NDSolve`ProcessEquations[Sequence @@ eq, options];
  count = 0;
  ctime = state@"CurrentTime"["Forward"];
  timer = First@Timing@While[ctime < tmax,
      Quiet[NDSolve`Iterate[state, tmax*1.000001], NDSolve::mxst];
      ctime = state@"CurrentTime"["Forward"];
      ic = state@"SolutionVector"["Forward"];
      dt = state@"TimeStep"["Forward"];
      state = First@NDSolve`Reinitialize[state, {x[ctime] == ic[[1]], x'[ctime] == ic[[2]]}, StartingStepSize -> dt];
      count++;
     ];
  {timer, state@"CurrentTime"["Forward"], count, state@"SolutionVector"["Forward"]}
]

(* ==> {3.1722, 1., 2948, {0.841471, 0.540302}} *)`

Timing results are not optimistic. The situation is worse when one want use step size from the previous step (with StartingStepSize -> dt commented in Reinitialize).

Is there a better (more efficient) way to to this?

share|improve this question
    
StepMonitor doesn't allow to change the current state of integrated equation. –  mmal Sep 28 '12 at 17:58
add comment

1 Answer

Not sure if this does what you want, but you can use StepMonitor or EvaluationMonitor and the trick to make it stop is to do a Throw[] first thing. You can obtain the state of the solution in there also

Catch[NDSolve[{y'[t] == 0, y[0] == 0}, y[t], {t, 0, 1}, 
    StepMonitor :> {Print[Row[{"t=", t, " y[t]=", y[t]}]]; Throw["one step done"]}]]

this prints the time of the one step and the value of y is then

t=0.000102648 y[t]=0.
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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