# Can't use my own discontinuous coefficient within NDSolve

I want to evaluate the pendulum's motion with changing equilibrium state given by eqilValue[t]:

tauArr = {50, 100, 150, 200, 250}; (* points of changing, there may be a lot of them *)
(* that's why i need a function *)
eqilValue[t_] := (If[0 <= t < tauArr[[1]], eV = 0];
If[t >= tauArr[[Length[tauArr]]], eV = -5];
For[j = 1, j <= (Length[tauArr] - 1), j++,
{If[tauArr[[j]] <= t < tauArr[[j + 1]], eV = j]}];
Return[eV]
);(* equilibrium states, for example *)

gamma = 0.1;
kappa = 0.1;

sol = First[
f /. NDSolve[{f''[t] + gamma*f'[t] + kappa*(f[t] - eqilValue[t]) == 0,
f[0] == 0, f'[0] == 0}, f, {t, 0, 300}]];
Plot[eqilValue[t], {t, 0, 400}]
Plot[sol[t], {t, 0, 300}, PlotRange -> All]


This code throws an error when launched first time, otherwise shows a wrong behavior, just like there is the only final state:

(*NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.. >>

ReplaceAll::reps: {NDSolve[{0.01 (Times[<<2>>]+f[<<1>>])+0.1 (f^\[Prime])[t]+(f^          \[Prime]\[Prime])[t]==0,f[0]==0,(f^\[Prime])[0]==0},f,{t,0,300}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>
*)


What am I doing wrong? Is there a better way to set up differential equations with discontinuities in the coefficients?

• Define eqilValue this way: eqilValue[t_?NumericQ] := (...). The need for and use of ?NumericQ is explained in this answer to the "Pitfalls" question and its links. – Michael E2 Nov 17 '14 at 21:15
• FWIW, you might be interested in Alternatives to procedural loops. Mathematica has tools built in to support Fortran/C/Java type programming, such as For, which lets folks transition more easily. But they don't work as well with the rest of M as other tools. – Michael E2 Nov 18 '14 at 13:00
• I took the liberty to edit the question to include the part about handling discontinuous function. If you don't like it, feel free to roll back the edit. – Michael E2 Nov 18 '14 at 14:13
• I'm voting to reopen. I think belisarius's answer shows there is an interesting aspect to the problem that the edited question now reflects. It also a good question for exhibiting some new abilities of DSolve. – Michael E2 Nov 18 '14 at 15:01

I believe it's better to cope with discontinuous functions in NDSolse[] by using WhenEvent[] and DiscreteVariables ->. In your case:

tauArr = {50, 100, 150, 200, 250};
whens = WhenEvent @@@ MapIndexed[{t == #1, ep[t] -> If[#2 == {5}, -5, #2[[1]]]} &, tauArr]
sol = NDSolve[{f''[t] + gamma*f'[t] + kappa*(f[t] - ep[t]) == 0,
whens, f[0] == 0, f'[0] == 0, ep[0] == 0}, f, {t, 0, 300},
DiscreteVariables -> {ep}]
Plot[f[t] /. sol[[1]], {t, 0, 300}, PlotRange -> All]


• This is a nice answer to a better question. I think I'll edit the question to match. – Michael E2 Nov 18 '14 at 12:45
• @MichaelE2 Thanks! :) – Dr. belisarius Nov 18 '14 at 14:03

NDSolve and now, in V10, DSolve, too, can handle differential equations with discontinuous coefficients. The coefficients have to be in terms of functions that will be recognized as discontinuous. Procedural programs such as in the definition of eqilValue above will not be recongnized and will be treated as numerical black boxes. Piecewise will be recognized and both NDSolve and DSolve analyze the function and attempt to determine where the discontinuities occur. One thing that is a little different with DSolve is that one has to ((always?) specify a finite interval. DSolve will construct a piecewise solution that is valid over the interval.

We can use either NDSolve or DSolve. I'll show DSolve. It is new and can actually give an "exact" solution -- exact, except for the approximate parameters gamma and kappa.

eqVal = Piecewise[
MapIndexed[{First[#2] - 1, t < #1} &, tauArr],
-5];

sol = DSolve[{f''[t] + gamma*f'[t] + kappa*(f[t] - eqVal) == 0, f[0] == 0,
f'[0] == 0}, f, {t, 0, 300}]


Plot[f[t] /. sol, {t, 0, 300}, PlotRange -> All]


Remarks: (1) If we replace DSolve by NDSolve above, we will of course get an interpolating function that closely approximates the exact solution. I'm not sure of the differences in how NDSolve processes the piecewise coefficient and how it processes the WhenEvent method of belisarius. Usually, NDSolve uses root-finding algorithms to locate events, but for equations of the form t == t0, root-finding and solving exactly are pretty much the same thing. They certainly are in this case where the t0 are integers. In any case, there's not a great difference in the two methods.

(2) The eqilValue function had a regular pattern that facilitated the construction of the Piecewise expression (and WhenEvent` expressions). For an irregular pattern, one can simply construct them by hand.

• Nice, didn't know about DSolve new juggles (Still on 9 here) – Dr. belisarius Nov 18 '14 at 14:08
• @belisarius Thanks. DSolve can also do delay differential equations (over a finite interval) in the same way (piecewise solution). – Michael E2 Nov 18 '14 at 14:11