7

Backslide introduced in 9.0, persisting through 11.0.1

I try to solve a boundary value problem of coupled ODEs. For reference, when I solve this set of equations it works without any problems:

eq1 = {y'[x] == Piecewise[{{z[x], x < 5}, {y[x], x > 5}}]};
eq2 = {z'[x] == Piecewise[{{2 z[x] + 1, x < 5}, {2 y[x] + 1, x > 5}}]};
eqs = {eq1, eq2, y[0] == 1, z[10] == 0};
sol1 = NDSolve[eqs, {y, z}, {x, 0, 10}, 
  Method -> {"Shooting", 
    "StartingInitialConditions" -> {y[0] == 1, z[0] == 1}}]

When I try to use more than 2 sections in the piecewise function I receive the following error:

eq1 = {y'[x] == 
    Piecewise[{{z[x], x < 5}, {y[x], 8 > x > 5}, {2 z[x], x > 8}}]};
eq2 = {z'[x] == Piecewise[{{2 z[x] + 1, x < 5}, {2 y[x] + 1, x > 5}}]};
eqs = {eq1, eq2, y[0] == 1, z[10] == 0};
sol1 = NDSolve[eqs, {y, z}, {x, 0, 10}, 
  Method -> {"Shooting", 
    "StartingInitialConditions" -> {y[0] == 1, z[0] == 1}}]

NDSolve::bvdisc: NDSolve is not currently able to solve boundary value problems with discrete variables. >>

I searched for a solution to this error online but non of the offered solutions helped. I don't understand why mathematica reads at as discrete variables when more than 2 sections are included. Same error appears whenever a boundary of the form a<x<b is set in the piecewise function. It is important to note that for initial boundary problems, there are no issues. For example:

eq1 = {y'[x] == 
    Piecewise[{{z[x], x < 5}, {y[x], 8 > x > 5}, {2 z[x], x > 8}}]};
eq2 = {z'[x] == Piecewise[{{2 z[x] + 1, x < 5}, {2 y[x] + 1, x > 5}}]};
eqs = {eq1, eq2, y[0] == 1, z[0] == 0};
sol1 = NDSolve[eqs, {y, z}, {x, 0, 10}]

Yields

{{y -> InterpolatingFunction[{{0., 10.}}, <>], 
  z -> InterpolatingFunction[{{0., 10.}}, <>]}}

Does anyone have a solution to this?

Thanks in advance.

xzczd
  • 65,995
  • 9
  • 163
  • 468
Eng.Technion
  • 73
  • 3

2 Answers2

9

It's a backslide introduced since v9:

enter image description here

enter image description here

and another issue of "DiscontinuityProcessing" so we at least have 2 workarounds.

One is to turn off the "DiscontinuityProcessing":

eq1 = {y'[x] == Piecewise[{{z[x], x < 5}, {y[x], 8 > x > 5}, {2 z[x], x > 8}}]};
eq2 = {z'[x] == Piecewise[{{2 z[x] + 1, x < 5}, {2 y[x] + 1, x > 5}}]};
eqs = {eq1, eq2, y[0] == 1, z[10] == 0};
sol1 = NDSolve[eqs, {y, z}, {x, 0, 10}, Method -> {"DiscontinuityProcessing" -> False}]

But as pointed out by Michael E2 in his answer, this is dangerous.

Another workaround is to transform the Piecewise into UnitStep, with the help of the undocumented Simplify`PWToUnitStep:

eq1 = {y'[x] == Piecewise[{{z[x], x < 5}, {y[x], 8 > x > 5}, {2 z[x], x > 8}}]};
eq2 = {z'[x] == Piecewise[{{2 z[x] + 1, x < 5}, {2 y[x] + 1, x > 5}}]};
eqs = {eq1, eq2, y[0] == 1, z[10] == 0};
sol1 = NDSolve[eqs // Simplify`PWToUnitStep, {y, z}, {x, 0, 10}]
Community
  • 1
xzczd
  • 65,995
  • 9
  • 163
  • 468
  • NDSolve implements discontinuities with discrete variables for Piecewise with more than two pieces, not counting the default. This approach may solve some problems, but causes one here. Why do you consider this a bug, since there is a proper warning, instead of it being a limitation? – Michael E2 Dec 14 '16 at 13:49
  • @MichaelE2 After a second thought I modified the tag to "backslide", given NDSolve doesn't output a wrong answer in this case. BTW you mean NDSolve is trying to handle the discontinuity with DiscreteVariables in this case? (At least in v9 the warning bvdisc isn't documented at all. ) Maybe you can elaborate a little with an answer? – xzczd Dec 14 '16 at 14:31
  • I won't have time, at least for several days. Here's a closely related answer: http://mathematica.stackexchange.com/a/106193/4999 -- I get the same behavior in V9 as V11 (and the message is not documented in V11, either). – Michael E2 Dec 14 '16 at 14:31
  • You can see a little of what's going on with {state} = NDSolve`ProcessEquations[<ndsolve code>] and with state["NumericalFunction"] and state["Variables"]. Image: http://i.stack.imgur.com/OWNZ8.png – Michael E2 Dec 14 '16 at 14:37
4

Interesting: Changing the conditions in the Piecewise function to simple one-sided inequalities allows NDSolve to work.

eq1 = {y'[x] == Piecewise[{{z[x], x < 5}, {2 z[x], x > 8}, {y[x], x > 5}}]};
eq2 = {z'[x] == Piecewise[{{2 z[x] + 1, x < 5}, {2 y[x] + 1, x > 5}}]};
eqs = {eq1, eq2, y[0] == 1, z[10] == 0};
sol1 = NDSolve[eqs, {y, z}, {x, 0, 10}, 
  Method -> {"Shooting", 
    "StartingInitialConditions" -> {y[0] == 1, z[0] == 1}}]

You still get a NDSolve::berr message (scaled boundary error), but that it not surprising in a discontinuous BVP.

Michael E2
  • 235,386
  • 17
  • 334
  • 747