NDSolve: ProcessEquations and Reinitialize with Piecewise functions

I am having trouble with using NDSolveReinitialize 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 @ NDSolveProcessEquations[simplesys, r, t]


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

newstate = NDSolveReinitialize[state, {r[0] == 2}]


this fails with the message:

NDSolveReinitialize::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?

• In this specific case you can replace Piecewise[...] with UnitStep[10 - t]. Aug 14, 2013 at 6:51
• Thanks for the quick response @b.gatessucks. Works absolutely fine now :) Aug 14, 2013 at 6:55

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@ NDSolveProcessEquations[simplesys, r, t,
Method -> {"EquationSimplification" -> "Residual"}];

newstate = First@NDSolveReinitialize[state, {r[0] == 2}];
NDSolveIterate[newstate, 20];
sol = NDSolveProcessSolutions[newstate];

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


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"

The OP's original code works as is in V10.0.0 through V11:

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

{state} = NDSolveProcessEquations[simplesys, r, t];


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

Module[{newstate, sol},
GraphicsRow@Table[
{newstate} = NDSolveReinitialize[state, {r[0] == r0}];
NDSolveIterate[newstate, 20];
sol = NDSolveProcessSolutions[newstate];
Plot[Evaluate[r[t] /. sol], {t, 0, 20}, AxesOrigin -> {0, 0}, PlotRange -> {0, 20}],
{r0, {0, 4, 6}}]
]


For V9.0.1, one may use

{state} = NDSolveProcessEquations[SimplifyPWToUnitStep@simplesys, r, t];


This is effectively the same thing as suggested in @b.gatessucks' comment.