# NDSolve with piece-wise function and BVP

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?

• 1. A ) is missing somewhere in equadiff, 2. What's rt? Please double check your code. Also, it's better to show us a complete sample reproducing the issue, rather than embed it in the text. Aug 18, 2022 at 2:10
• Thanks, edited. Aug 18, 2022 at 9:28

(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\}$$

• Thanks, edited away parameters (not useful here). I want to solve numerically because I want to solve a more general case that is not solvable analytically. Aug 17, 2022 at 20:52
• That solution doesn't appear to satisfy m[100]==0. Aug 18, 2022 at 6:56
• @BillWatts Why do you say it doesn't satisfy m[100]==0? Aug 18, 2022 at 19:12
• Plug in 100 for t. It's small but not zero. Aug 18, 2022 at 21:23
• No, look at the equation. It is a standard exponential with negative exponent that asymtotically approaches zero at infinity, not 100, even if you use 1/2. Aug 18, 2022 at 23:33

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 SimplifyPWToUnitStep instead of "DiscontinuityProcessing" -> False:

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