How to use NDSolve to track a deviation from equilibrium?

Question

i have a task to find a event, to achieve this i think need to integrate until a constant function suddenly breaks the equilibrium and give a sharp peak. I think that EventLocator should work, thus tried the proposals from How to use NDSolve to track equilibrium?, with a threshold, but the numerical calculator is still complaining that the value before and after the step are not real numbers, but this would not be the case in theory, the code:

psi=Pi;
For[i=1,i<141,i++,w[i]=1/2];
fext=1;
k = 35.529;
ci = Table[y[i][0] == RandomReal[{-Pi, Pi}], {i, 1, 141}];
AppendTo[ci,y[0][0] == y[141][0] - psi];

lckm = {y[0]'[t] == fext, y[141]'[t] == fext,Table[y[i]'[t] == 
         w[i] + 
          k (Sin[y[i + 1][t] - y[i][t]] + 
             Sin[y[i - 1][t] - y[i][t]]), {i, 1, 
         140}]};

sol = NDSolve[{lckm, ci}, 
  Table[y[i], {i, 0, 141}], {t, 0.0, 800.0},
  WorkingPrecision -> 30, AccuracyGoal -> 15, PrecisionGoal -> 15, 
  MaxSteps -> Infinity];

v1 = y[140][t] /. sol[[1]];
v2 = y[141][t] /. sol[[1]];
norm = Sin[v1 - v2] /. t -> 800;
k = 35.526;
    ci = Table[
       y[i][0] == 
        Part[{y[i][t] /. sol} /. t -> 800.0, 1], {i, 1, 141}];
    AppendTo[ci, y[0][0] == y[141][0] - psi];


lckm = {y[0]'[t] == fext, y[141]'[t] == fext,Table[y[i]'[t] == 
     w[i] + 
      k (Sin[y[i + 1][t] - y[i][t]] + 
         Sin[y[i - 1][t] - y[i][t]]),{i, 1, 
     140}]};
sol = NDSolve[{lckm, ci}, 
  Table[y[i], {i, 0, 141}], {t, 0.0, 800.0}, 
  Method -> {"EventLocator", 
    "Event" -> 
     Boole[Abs[
        Sin[y[140][t] - y[141][t]]/norm - 1] > 
       1/10], "EventAction" -> Print["t:", t]}, 
  WorkingPrecision -> 30, AccuracyGoal -> 15, PrecisionGoal -> 15, 
  MaxSteps -> Infinity]

All the constant have a value well ascertained value, including norm that is a normalization responsable for the deviation detect. The answer to this, in mathematica 7, is:

NDSolve::precw: The precision of the differential equation ({((y[0])^\[Prime])[t]==1,((y[141])^\[Prime])[t]==1,((y[1])^\[Prime])[t]==1/2+66.439 (Sin[Subscript[\[Null], <<2>>][<<1>>]+Times[<<2>>]]-Sin[Plus[<<2>>]]),<<6>>,((y[8])^\[Prime])[t]==1/2+66.439 (Sin[\[Null], <<2>>][<<1>>]+Times[<<2>>]]-Sin[Plus[<<2>>]]),<<274>>}) is less than WorkingPrecision (30.`). >>

NDSolve::evre: The value of the event function at t = 2.79284074034703538067228005241`30.*^-10 was not a real number.  The event will be ignored in steps where it does not evaluate to real numbers at both ends. >>

Please note that this logic is set to work, i.e. this error is not dependent in the value of norm in the sense that it needs it to detect the deviation. The point is: why Mathematica is outputing an error saying that something is not a real number ? where is the thing that is not a real number ?

Iam sorry if the example is too extense, i think that the problem poses it in this form or iam too newbie in mathematica and dont know rewrite.

Does anyone has any suggestion in order to circumvent it ?

Thank you !!!

Answer 1

Now, you did put some pressure on me by giving a working example :-). This is most probably not a solution yet, but some thoughts after trying to run your code:

Problem Size: that's another point in making your question more friendly to people trying to answer it: make the problem size as small as possible, it's just no fun to wait minutes for a result just to see it's structure. I reduced your example from 141 equations to 11, which will of course be much simpler to handle. In this case the 141 equations are solved very fast with machine precision, but it took too long to run these with the high precision that you define, so I stopped it and reduced the problem size. Trying new things as your event handler to detect deviation from equilibrum is also something I'd strongly recommend to first try with a problem as simple as possible. Only when you have that working alright you should run your real problem.

WorkingPrecision: since you give it equations with lower precision than what you want NDSolve to use it complains (thats the NDSolve::precw warning). What it basically does (I think) is giving that warning and ignore it. If you want it to really use that precision, you'll need to give your input in at least that precision. What I tried was making the input exact with Rationalize[...,0] which indeed will not give the NDSolve::precw warning anymore. I don't think you'll actually need it, at least for my reduced example the (trivial) result doesn't change when using MachinePrecision but of course will be much faster...

Real Number: the warning about something not being a real number is issued because the right hand side of the initial conditions ci is not set to a number but to a list with one number, which you can see when looking e.g. at ci[[3]]. It's easy enough to change that by adding another 1 to Part. With these changes the code will work in Version 8, but not in Version 7. I think this is probably because of an error which will replace the integers in y[n][t] to real numbers, so in the event checks there appears something like y[20.][0.0144971] which doesn't work. One way to circumvent that is to replace those integers with strings, like I show in the code below. For other readers I'd like to emphasize that this is only necessary for Version 7 (and maybe earlier versions as well).

Version Information: I admit that in this case this was absolutely not to be expected, but part of your problem seem to be a bug in version 7. Noone trying with version 8 could have ever see that error, so this information was essential to have a chance to help you, so if not using the most recent version it probably is always useful to include that information as you did.

RuleDelayed for event specification: You will notice that I have used RuleDelayed (shortcut: :>) instead of Rule (shortcut ->) to define the values for the option "EventAction" as well as "Event". This is probably only necessary for "EventAction" always, but usualy is also a good idea for "Event": you want these expression only to be evaluated at every single step and the RuleDelayed will ensure that this is so.

Altogether the following will work with no errors (but doesn't stop because the problem doesn't actually show the behavior that the event tries to tests for):

psi = Pi;
Do[w[i] = 1/2, {i, 1, 11}];
fext = 1;
k = 35529/1000.;
ci = Table[y[i][0] == RandomReal[{-Pi, Pi}], {i, 1, 11}];
AppendTo[ci, y[0][0] == y[11][0] - psi];

lckm = {y[0]'[t] == fext, y[11]'[t] == fext, 
   Table[y[i]'[t] == 
     w[i] + k (Sin[y[i + 1][t] - y[i][t]] + 
         Sin[y[i - 1][t] - y[i][t]]), {i, 1, 10}]};

sol = NDSolve[{lckm, ci}, Table[y[i], {i, 0, 11}], {t, 0, 800}];

v1 = y[10][t] /. sol[[1]];
v2 = y[11][t] /. sol[[1]];
norm = Sin[v1 - v2] /. t -> 800;
k = 35526/1000.;
ci = Table[
   y[i][0] == Part[{y[i][t] /. sol} /. t -> 800, 1, 1], {i, 1, 11}];
AppendTo[ci, y[0][0] == y[11][0] - psi];

lckm = {
   y[0]'[t] == fext,
   y[11]'[t] == fext, 
   Table[y[i]'[t] == 
     w[i] + k (Sin[y[i + 1][t] - y[i][t]] + 
         Sin[y[i - 1][t] - y[i][t]]), {i, 1, 10}]
   };

sol = NDSolve[Evaluate[{lckm, ci} /. y[n_] :> y[ToString[n]]], 
   Table[y[ToString[i]], {i, 0, 11}], {t, 0.0, 800}, 
   Method -> {"EventLocator", 
     "Event" :> (Abs[Sin[y["10"][t] - y["11"][t]]/norm - 1] - 1/10
       ), "EventAction" :> Print["t:", t]}, MaxSteps -> Infinity];

Plot[Evaluate[Table[y[ToString[k]][t] /. sol, {k, 1, 11}]], {t, 0, 800}]

If your real problem needs it and you want NDSolve to really use WorkingPrecision->30, then you could run the NDSolve as here:

sol = NDSolve[Evaluate[Rationalize[{lckm, ci}, 0]], 
   Table[y[i], {i, 0, 11}], {t, 0, 800},
   WorkingPrecision -> 30, AccuracyGoal -> 15, PrecisionGoal -> 15, 
   MaxSteps -> Infinity
   ];

which will not issue any warnings and should really use the higher precision, but of course will take a lot longer than the standard settings (I think it will switch to use software floating point arithmetic for these settings, which is a lot more expensive...).

---

PS: (you'll note that I coudn't resist to get rid of the For-loop, for reasons I gave in Are there any cases when For[] loops are reasonable?).