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 am having trouble with using NDSolve`Reinitialize when the system consists of a piecewise function.

If we define the ODE system

simplesys = {r'[t] == Piecewise[{{1, 0<= t <=10}, {0, 10<= t <=20}}, 0], r[0] == 0};

and process the equations with

state = First @ NDSolve`ProcessEquations[simplesys, r, t]

this works perfectly fine. However, when I try to reinitialize the ODE system with a new initial condition:

newstate = NDSolve`Reinitialize[state, {r[0] == 2}]

this fails with the message:

NDSolve`Reinitialize::ntcs: Cannot solve constraint equations for initial conditions.

This is due to the piecewise function in the system. I don't have this problem for other types of systems, and it only occurs when I have a piecewise function. How can I solve this issue?

share|improve this question
3  
In this specific case you can replace Piecewise[...] with UnitStep[10 - t]. –  b.gatessucks Aug 14 '13 at 6:51
    
Thanks for the quick response @b.gatessucks. Works absolutely fine now :) –  jaclea Aug 14 '13 at 6:55

2 Answers 2

It is quite tricky! Piecewise[ ] functions work only with "EquationSimplification" -> "Residual"... I'll try to dig up why

simplesys = {r'[t] == Piecewise[{{1, 0 <= t <= 10}, {0, 10 <= t <= 20}}, 0], r[0] == 0};
state = First@ NDSolve`ProcessEquations[simplesys, r, t, 
                                  Method -> {"EquationSimplification" -> "Residual"}];

newstate = First@NDSolve`Reinitialize[state, {r[0] == 2}];
NDSolve`Iterate[newstate, 20];
sol = NDSolve`ProcessSolutions[newstate];

Plot[Evaluate[r[t] /. sol], {t, 0, 20}, AxesOrigin -> {0, 0}]

Mathematica graphics

The truth is that the "Residual" method doesn't seem to do much work, as it isn't simplifying the equations. I'm not sure if using the pre-processing scheme is a good idea if you're forced to use "Residual"

share|improve this answer

I think this might have been a bug. At least as of V10.0.0, the following works without incident:

simplesys = {r'[t] == Piecewise[{{1, 0 <= t <= 10}, {0, 10 <= t <= 20}}, 0], r[0] == 0};

{state} = NDSolve`ProcessEquations[simplesys, r, t];

Here we reinitialize three times and plot the results of each:

Module[{newstate, sol},
 GraphicsRow@Table[
   {newstate} = NDSolve`Reinitialize[state, {r[0] == r0}];
   NDSolve`Iterate[newstate, 20];
   sol = NDSolve`ProcessSolutions[newstate];
   Plot[Evaluate[r[t] /. sol], {t, 0, 20}, AxesOrigin -> {0, 0}, PlotRange -> {0, 20}],
   {r0, {0, 4, 6}}]
 ]

Mathematica graphics

share|improve this answer

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.