NDSolve with piece-wise function and BVP

Question

How to numerically solve a system of differential equation with boundary conditions and piece wise affine functions ?

Consider the following system of differential equations from some optimal control problem:

δ = Piecewise[{{0.0105, 0 <= t <= 10}, {0.0413,10 <= t <= 80}, {0.001, 80 <= t}}] ; p=Exp[-δ*0.05*t];    
equadiff = {m'[t] == 0.5 m[t] - (v[t]/(p * δ)^(-1/2)), v'[t] == -0.5v[t]}

Running

solution = NDSolve[Join[equadiff, {m[0] == 100}, {m[100]==0}], {m[t], v[t]}, {t, 0, 100}]

I get the following error:

I also get the same error when using the shooting method with initial value on variable v

I tried suggestions in this post on discontinuous data

When running

solution = NDSolve[Join[equadiff, {m[0] == 100}, {m[100]==0}], {m[t], v[t]}, {t, 0, 100},Method ->{"PDEDiscretization" -> "FiniteElement"}]

I get the following error

When running

solution = NDSolve[Join[equadiff, {m[0] == 100}, {m[100]==0}], {m[t], v[t]}, {t, 0, 100},Method -> {"DiscontinuityProcessing" -> False}]

I get the following errors:

Using finite elements tries to solve the system but the errors suggest that it encounters problem from divided by 0 or negative roots. I suspect that variable v somehow is numerically negative or 0 at some time points (which should be theoretically strictly positive in my problem).

How to solve this system with boundary conditions and piecewise data? I want to solve numerically, not analytically because it is particular case of a more general problem with no analytical solution. What other methods could I try?

Answer 1

(You haven't specified what parameters are.)

In this case, DSolve gives an analytical solution.

DSolve[Join[{m'[t] == 0.5 m[t] - (v[t]/(p*delta)^(-1/2)), v'[t] == -0.5 v[t]}, 
    {m[0] == 100}, {m[100] == 0}], {m[t], v[t]}, {t, 0, 100}]//Chop

$$\left\{\left\{m(t)\to 100. e^{-0.5 t},v(t)\to \frac{100. e^{-0.5 t}}{\sqrt{delta\ p}}\right\}\right\}$$

Answer 2

You're almost there. Just add a high enough WorkingPrecision:

solution = 
 NDSolveValue[
  Rationalize[#, 0] &@{equadiff, m[0] == 100, m[100] == 0}, {m, v}, {t, 0, 100}, 
  Method -> {"DiscontinuityProcessing" -> False}, WorkingPrecision -> 48]

ListLinePlot[solution, PlotRange -> All]

I've used an documented syntax of ListLinePlot here, see this post for more info.

BTW, you can also use Simplify`PWToUnitStep instead of "DiscontinuityProcessing" -> False:

solution = 
 NDSolveValue[
  Rationalize[#, 0] &@{equadiff, m[0] == 100, m[100] == 0} // 
   Simplify`PWToUnitStep, {m, v}, {t, 0, 100}, WorkingPrecision -> 48]