NDSolve problems of convergence

Question

I have a big function that I have to maximize, so I have to evaluate some equation with thousands of different values of the parameter. But the FindMaximum returns me many errors. I have traced the cause down to a differential equation that is part of the calculation. This differential equation has to be solved as function of the parameter (called "sup" in this case) and even though it is a seemingly "easy" function and I don't need a very high precision, in many points it does not converge and reaches the maximum number of steps. The following is a simplified sample code to view the problem with a Manipulate: as you slide the cursor with the Sup variable the NDsolve produces errors.

     Manipulate[
 Module[{pend, tmax}, tmax = 30; 
  pend = NDSolve[{1.1089 α''[t] + 1.2936 Sin[α[t]] + 
       210*sup*
        (0.33 Sin[α[t]] + (
          0.33^2 Cos[α[t]] Sin[α[t]])/Sqrt[
          1 - 0.33^2 Sin[α[t]]^2]) *
        UnitStep[Sin[α[t] + Pi*UnitStep[-α'[t]]]] == 
      0, α[0] == -Pi/3, α'[0] == 0}, α, {t, 0, 
     tmax}, MaxSteps -> 10^5, PrecisionGoal -> 8, 
    AccuracyGoal -> acc]; 
  Plot[α[x] /. pend, {x, 0, tmax}]], {{sup, 0.01, "Sup"}, 0.01,
   0.1, Appearance -> "Labeled"}, {{acc, 8, "Acc"}, 6, 10, 1, 
  Appearance -> "Labeled"}, Bookmarks -> {
   "Error1" :> {sup = 0.0577, acc = 7}, 
   "Error2" :> {sup = 0.0413, acc = 8}, 
   "Error3" :> {sup = 0.0351, acc = 6}}]

I have tried with increasing the MaxSteps and reducing the AccuracyGoal/PrecisionGoal, but it seems that there are always points in which the NDSolve does not return a good solution. And thus I cannot execute the FindMaximum to maximize my function.

Any hint on how to avoid this kind of problem and get NDSolve to return a correct result? Is there a way to trace down why NDSolve has this problems of convergence with this function?

Thanks for any help,

JBB

PD: Code update to define tmax. In the manipulate Bookmarks I have placed some examples of errors.

Answer 1

I finally found a "dirty" patch that I'm using and I leave here in case it is of use to someone else. Augmenting the WorkingPrecision is a good solution but makes the program extremely slow, thanks to xzczd in any case.

So I'm using some "knowledge" of the function:

• I know that the function is a decreasing oscillation and if it "shoots up" it is because of an integration problem, but not a real data,
• and also as I don't need a very high precision

What I do is to detect that there has been an integration problem (because the integration time does not reach the specified one) and if this is the case, then change a little bit the input parameter and repeat the integration until I get a "good" solution. It's not an elegant solution, but it works for me.

Here is the code:

    Module[{pend, tmax, supTemp},
 Manipulate[
  supTemp = sup;
  Label[repeatChangingSup]; (*If NDSolve has errors returns here and \
tries again with a slightly different surface *)
  tmax = 100;
  pend = Quiet@
    NDSolve[{1.1089 α''[t] + 1.2936 Sin[α[t]] + 
        210*supTemp*
         (0.33 Sin[α[t]] + (
           0.33^2 Cos[α[t]] Sin[α[t]])/Sqrt[
           1 - 0.33^2 Sin[α[t]]^2]) *
         UnitStep[Sin[α[t] + Pi*UnitStep[-α'[t]]]] == 
       0, α[0] == -Pi/3, α'[0] == 0}, α, {t, 0, 
      tmax}];
  If[tmax == pend[[1]][[1]][[2]][[1]][[1]][[2]],(* 
   If tmax is the same as the last value of the interp function, 
   the NDSolve has finished succesfully*)
   pend, (*if everything OK return pend*)
   supTemp = supTemp*1.0001; Print[supTemp]; 
   Goto[repeatChangingSup];]; (* 
  Else there was an error in NDSolve and we try with a sligthly \
different sup *)
  Plot[α[x] /. pend, {x, 0, tmax}],
  {{sup, 0.01, "Sup"}, 0.01, 0.1, Appearance -> "Labeled"}]]

Answer 2

A higher WorkingPrecison will somewhat help:

Manipulate[Module[{pend, tmax}, tmax = 30;
  pend = Quiet@
    NDSolve[{1.1089 α''[t] + 1.2936 Sin[α[t]] + 
        210*sup*(0.33 Sin[α[
              t]] + (0.33^2 Cos[α[t]] Sin[α[t]])/
            Sqrt[1 - 0.33^2 Sin[α[t]]^2])*
         UnitStep[Sin[α[t] + Pi*UnitStep[-α'[t]]]] == 
       0, α[0] == -Pi/3, α'[0] == 0}, α, {t, 0, 
      tmax}, MaxSteps -> 10^5, PrecisionGoal -> 8, 
     AccuracyGoal -> acc, WorkingPrecision -> 32];
  Plot[α[x] /. pend, {x, 0, tmax}]], {{sup, 0.01, "Sup"}, 0.01,
   0.1, Appearance -> "Labeled"}, {{acc, 8, "Acc"}, 6, 10, 1, 
  Appearance -> "Labeled"}, 
 Bookmarks -> {"Error1" :> {sup = 0.0577, acc = 7}, 
   "Error2" :> {sup = 0.0413, acc = 8}, 
   "Error3" :> {sup = 0.0351, acc = 6}}]

Notice that here I use Quiet to quit the warning message because the warning generates just for those decimates (their Precisions are MachinePrecision) in the equation and it doesn't influence the result. If the equation isn't in Manipulate you can also Rationalize all the decimates in the equation to avoid the warning.

For more information about the reason why I change WorkingPrecision, you can have a look at the following questions (there may be more, but at the moment I just remember these):

ReplaceAll[] and Limit[] don't give correct results for this expression under extreme variables

Funny behaviour when plotting a polynomial of high degree and large coefficients

FindRoot errors

NSolve gives additional solutions that don't satisfy the equations!