Can't use my own discontinuous coefficient within NDSolve

Question

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?

Answer 1

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]

Answer 2

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.