1
$\begingroup$

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:

NDSolve 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

enter image description here

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:

enter image description here enter image description here

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?

$\endgroup$
2
  • 1
    $\begingroup$ 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. $\endgroup$
    – xzczd
    Aug 18 at 2:10
  • $\begingroup$ Thanks, edited. $\endgroup$
    – Kredan
    Aug 18 at 9:28

2 Answers 2

1
$\begingroup$

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

$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$
    – Kredan
    Aug 17 at 20:52
  • 1
    $\begingroup$ That solution doesn't appear to satisfy m[100]==0. $\endgroup$
    – Bill Watts
    Aug 18 at 6:56
  • $\begingroup$ @BillWatts Why do you say it doesn't satisfy m[100]==0? $\endgroup$ Aug 18 at 19:12
  • 1
    $\begingroup$ Plug in 100 for t. It's small but not zero. $\endgroup$
    – Bill Watts
    Aug 18 at 21:23
  • 1
    $\begingroup$ 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. $\endgroup$
    – Bill Watts
    Aug 18 at 23:33
1
$\begingroup$

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]

enter image description here

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]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.