3
$\begingroup$

Good afternoon,

I am attempting to solve the following ODE:

d[h_] := Piecewise[{{0, Abs[h] > 110}, {0.1, 
110 >= Abs[h] > 100 && h < 0}, {-0.1,  
110 >= Abs[h] > 100 && h >= 0}}]

S = NDSolve[{h''[t] == d[h[t]], h'[-120] == 1, h[-120] == -120}, h, {t, -120, 120}]

From the definition of d[h_], we would expect h'[t] to be constant in the regions where d[h_] = 0, and otherwise linearly increasing/decreasing depending on whether h[t] is positive or negative. And indeed, using the bounds -120 to 120, we see that that is the case:

Plot[h[t] /. S, {t, -120, 120}, AxesLabel -> {"t", "h[t]"}, PlotLabel -> "tf=120"]

h[t]

Plot[h'[t] /. S, {t, -120, 120}, AxesLabel -> {"t", "h'[t]"}, PlotLabel -> "tf=120"]

h'[t] From my knowledge of the analytical solution, this appears to be the correct answer.

Now, let me take tf=120->300, but examine h[t] and h'[t] only in the region t=-120 to 120, the same as before. The results dramatically change:

S2 = NDSolve[{h''[t] == d[h[t]], h'[-120] == 1, h[-120] == -120}, h, {t, -120, 300}]

Plot[h[t] /. S2, {t, -120, 120}, AxesLabel -> {"t", "h[t]"}, PlotLabel -> "tf=300"]

h[t] with tf=300

Plot[h'[t] /. S2, {t, -120, 120}, PlotRange -> All, 

AxesLabel -> {"t", "h'[t]"}, PlotLabel -> "tf=300"]

h'[t] with tf=300

The answer is now incorrect!

It appears that NDSolve is somehow not changing h'[t] when h''[t] is non-zero (for t such that 100<=h[t]<110). However, I'm confused why this is the case when all I did was change tf.

My question is that is there a way to prevent this from occurring? This is a simplified version of a much more complicated problem I am working on, which does the same thing.

I have tried making the following change to no avail:

S2 = NDSolve[{h''[t] == d[h[t]], h'[-120] == 1, h[-120] == -120}, h, {t, -120, 300}, InterpolationOrder -> All, WorkingPrecision -> 22]

and have also played around with the option:

S2 = NDSolve[{h''[t] == d[h[t]], h'[-120] == 1, h[-120] == -120}, h, {t, -120, 300}, InterpolationOrder -> All, WorkingPrecision -> 22, Method -> {"DiscontinuityProcessing" -> False}]

but neither change the last two plots at all. I'm not sure what else to try/any potential workarounds.

Thanks!

$\endgroup$
5
$\begingroup$

NDSolve gets excited when the estimated error is zero, which happens when solutions are of low degree. Rein it in with MaxStepSize or MaxStepFraction:

S2 = NDSolve[{h''[t] == d[h[t]], h'[-120] == 1, h[-120] == -120}, 
  h, {t, -120, 300}, MaxStepSize -> 1]

Plot[h[t] /. S2, {t, -120, 120}, AxesLabel -> {"t", "h[t]"}, 
 PlotLabel -> "tf=300"]

enter image description here

|improve this answer|||||
$\endgroup$
  • 2
    $\begingroup$ (+1) just for NDSolve gets excited ;) $\endgroup$ – Nasser Mar 25 at 23:59
  • $\begingroup$ @Nasser Thanks, I guess I got excited and left the sentence incomplete. :) $\endgroup$ – Michael E2 Mar 26 at 0:01
  • $\begingroup$ Thank you so much, that fixed my problem both on the simplified and full version of my problem! $\endgroup$ – leafyrabbit Mar 26 at 18:12

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.