1
$\begingroup$

I'm trying to solve a set of ODEs with a Heaviside step function which depends on the sign of the derivative of the function.

This is a simplifying example of what I'm trying to do

sol = NDSolve[{
x'[t] == (0.1 x[t]^3 - x[t]) HeavisideTheta[x'[t]], 
x[0] == -2}, {x}, {t, 0, 1}];
Plot[Evaluate[x[t]] /. sol, {t, 0, 10}]

I get a lot of errors (NDSolve::tddisc, NDSolve::ntdvdae, NDSolve::nlnum), and it seems that it is not the proper way to work with NDSolve.

Can anyone please advise on how to properly use the Heaviside function inside NDSolve?

|improve this question
$\endgroup$
  • 2
    $\begingroup$ Use UnitStep instead. $\endgroup$ – AccidentalFourierTransform Oct 23 '18 at 19:30
  • 1
    $\begingroup$ Heaviside -> UnitStep does not help for me (V11.3). In both cases I get a solution -- is the solution correct? I think so. It's impossible to solve for x'[t], I think, which makes it a DAE, which NDSolve can solve. This equivalent: NDSolve[{x'[t] == Clip[0.1 x[t]^3 - x[t], {0, Infinity}], x[0] == -1/2}, {x}, {t, 0, 1}] but emits no warnings or errors. $\endgroup$ – Michael E2 Oct 24 '18 at 0:41
  • 1
    $\begingroup$ Those warnings may or may not mean your results are invalid. Your code does get results that that need checking, but you do need to use the same range on x in Plot as you use in NDSolve $\endgroup$ – Bill Watts Oct 24 '18 at 7:46
  • $\begingroup$ You can plot x'[t] $\endgroup$ – Bill Watts Oct 24 '18 at 7:48
2
$\begingroup$

I guess your main question is if your solution is correct. One way to do it:

sol = NDSolve[{x'[t] == (0.1 x[t]^3 - x[t]) HeavisideTheta[x'[t]], 
    x[0] == -2}, x[t], {t, 0, 10}] // Flatten;

And yes, I get the same warnings you do, although I get one less warning with UnitStep but I get the same answer.

x[t_] = x[t] /. sol;

Plot[Evaluate[x[t]], {t, 0, 10}, PlotRange -> All]

enter image description here

It satisfies the condition at x == 0. And to check if it satisfies the diffeq:

Plot[Evaluate[x'[t] - (0.1 x[t]^3 - x[t]) HeavisideTheta[x'[t]]], {t, 0, 10}, PlotRange -> All]

enter image description here

Looks reasonably close to me.

|improve this answer
$\endgroup$
1
$\begingroup$

As already mentioned by the comments and previous answer, the warning doesn't necessarily mean the result is unreliable. If you still feel uncertain, the following is another way to avoid the warning.

We first solve for $x'(0)$:

eq = x'[t] == (0.1 x[t]^3 - x[t]) HeavisideTheta[x'[t]];

newic = Equal @@@ 
  First@Solve[eq /. HeavisideTheta -> UnitStep /. t -> 0 /. x[0] -> -2, x'[0]]

(* {x'[0] == 1.2} *)

Then define an approximate HeavisideTheta:

appro = With[{k = 1000}, ArcTan[k #]/Pi + 1/2 &];

Replace the HeavisideTheta with appro and differentiate the equation once:

neweq = D[eq /. HeavisideTheta -> appro, t]

And you'll find NDSolve can solve the new equation without any warning:

sol = NDSolve[{neweq, x[0] == -2, newic}, x, {t, 0, 10}];
Plot[Evaluate[x[t]] /. sol, {t, 0, 10}, PlotRange -> All]

The resulting picture looks the same as that in Bill Watts's answer so I'd like to omit it here.

|improve this answer
$\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.