# Solving a two point boundary problem for a piecewise system of equations

I'm trying to solve a two point boundary value problem with a piecewise system of equations like:

dx[t_] := Piecewise[{{0, y[t] <= 0}, {x[t]^(3/4) - x[t], y[t] > 0}}]
dy[t_] := Piecewise[{{0, y[t] <= 0}, {(1 - y[t])/x[t]^(1/4) + y[t], y[t] > 0}}]

NDSolve[{x'[t] == dx[t], y'[t] == dy[t], x == 1, y == 2}, {x, y}, {t, 0, 10}]


which produces the error

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


I tried the answers here by entering

funs = ParametricNDSolveValue[{x'[t] == dx[t], y'[t] == dy[t], x == 1, y == yt}, {x, y}, {t, 0, end}, yt]

FindRoot[funs[yt][end], {yt, 2}]


which yields the same error.

Is there any way to solve this system? Thanks for the help.

There is no need to use Piecewise because both $x(t)$ as well as $y(t)$ remain non negative. Furthermore, the equation for x is independent of y and can be solved analytically.

The result $x(t)$ is then put into the equation for $y(t)$ which will be solved numerically.

Finally we verify our first statement about the non-negativity.

Ok, here we go.

The eqution for $x(t)$ is

sx = DSolve[x'[t] == x[t]^(3/4) - x[t] && x == 1, x[t], t]


We can safely ignore the error messages; they are just information which is not relevant here.

(* Out = {{x[t] -> E^-t (-2 + E^(t/4))^4}, {x[t] -> E^-t ((-1 - I) + E^(t/4))^4}, {x[t] -> E^-t ((-1 + I) + E^(t/4))^4}} *)


Taking the real-valued solution

xx[t_] = x[t] /. sx[]

(* Out= E^-t (-2 + E^(t/4))^4 *)

Plot[xx[t], {t, 0, 30}] Now solving numerically the equation for $y(t)$ using the previously obtained function $xx(t)$:

sy = NDSolve[y'[t] == (1 - y[t])/xx[t]^(1/4) + y[t] && y == 10, y[t], {t, 0, 50}][]


$\{y[t]\to \text{InterpolatingFunction}[\{\{0.,50.\}\},<>][t]\}$

Works out smoothly, no precision problems. Here's the function itself

yy[t_] = y[t] /. sy


$\text{InterpolatingFunction}[\{\{0.,50.\}\},<>][t]$

Finally ploting both functions together, magnifying $x(t)$ by a factor of ten for better visibility:

Plot[{10 xx[t], yy[t]}, {t, 0, 30}, PlotRange -> {-1, 20}, PlotLabel -> "blue x(t), red y(t)"] Nice problem. Hope this helps.

Regards, Wolfgang

• Thanks Wolfgang. That almost seemed to work. However, the two point boundary problem remains. Notice that x == 1, y == 2. With this, y(t) can be negative and thus the equation for x becomes dependent on y. – Mauricio Sep 9 '14 at 7:14
• @Mauricio: you are right. I took different initial conditions. Tried to reconsider the problem defining a time t0 below which y is zero. We then have x(t) = 1 for t<t0 and x(t) = E^(-t + t0)*(-2 + E^((t - t0)/4))^4 for t>t0. Assuming a value for t0, we can NDSolve for y. But I found no value of t0 which would lead to a consistent result. Mabe we can prove that no solution exists for the given initial conditions. – Dr. Wolfgang Hintze Sep 9 '14 at 8:54
• @Mauricio: I have found a consistent solution, see my 2nd answer. – Dr. Wolfgang Hintze Sep 9 '14 at 11:46

Here is now the - hopefully - correct (numeric) solution of the original problem.

To begin with, we define a time t0 such that y = 0 for 0<=t<=t0.

Für x(t) we have x(t) = 1 for 0<=t<=t0, otherwise we have first to solve the equation

sx = DSolve[x'[t] == x[t]^(3/4) - x[t] && x == 1, x[t], t]


We can ignore the error Messages.

(* Out= {{x[t] -> E^-t (-2 + E^(t/4))^4}, {x[t] ->
E^-t ((-1 - I) + E^(t/4))^4}, {x[t] -> E^-t ((-1 + I) + E^(t/4))^4}} *)


Taking the real-valued solution:

x[t] /. sx[]

(* Out= E^-t (-2 + E^(t/4))^4 *)


Having caluclated that, we now define the complete function x[t,t0]

x[t_, t0_] := E^(-t + t0) (-2 + E^((t - t0)/4))^4 /; t > t0

x[t_, t0_] := 1 /; 0 <= t <= t0


Plotting it (for t0 = 5)

Plot[{x[t, 5]}, {t, 0, 30}, PlotRange -> {0, 1}] Now inserting x[t,t0] we can NDSolve for y(t). Here we have selected the parameter t0 by trial and error such that y=2

t0 = 1.03214;
ta = 0;
te = 12;
sn = NDSolve[{D[y[t], t] == 1 - y[t] (1/x[t, t0]^(1/4) - 1), y[t0] == 0},
y[t], {t, ta, te}][];
{ta, te} = Head[(List @@ sn[])[]][[1, 1]];
yy[t_] = y[t] /. sn;
Print["y = ", yy];
Plot[{2, Evaluate[x[t, t0]], yy[t]}, {t, ta, te}, PlotRange -> {-1, 2.2},
PlotLabel ->
"Solution of the two-point boundary problem\nred curve = x(t), green curve \
= y(t)\nt0 = " <> ToString[t0] <> "\ty(10) = " <> ToString[yy]]

y = 2. Of course, y(t) must be set to zero below t=t0. This is trivial and not shown in the picture.

This seems to be the only solution. For greater values of t0 the boundary value y(10) decreases, and for smaller it increases.

Regards, Wolfgang

• Honestly speaking, I don't have much confidence in this solution. I don't like the sharp kink in y which comes from the zero in x ... – Dr. Wolfgang Hintze Sep 9 '14 at 12:28
• Thanks a lot Wolfgang. It does seems to work. I'll try to think if the kink is right or not. Cheers, Mauricio – Mauricio Sep 9 '14 at 13:11