Updated Lobatto

[HansOlsson]

Solution outlined in #4800
Might add a more detailed description as a comment.

The basic ideas are:

• Use tolerance in a somewhat normal way, instead of the "clever" is+delta==is trick.
  • Note: Originally there was a division by 10*Modelica.Constants.eps - the eps is removed but the 10 remains.
• Guard against oscillating integral being close to zero, by integrating the absolute of the function as well, and using that for comparison. Currently it uses max(abs(is), isabs/10) - we could change the factor 10 - or even just use isabs. It's a trade-off.

[HansOlsson]

As for the epsilion - originally it was

   is=x...*tol/(10*epsilon)
...
   if is+(i1-i2)==is

Not being clever means rewriting this as:

   is=x...*tol/(10*epsilon)
...
   if abs(i1-i2)<=abs(is)*epsilon

Which is basically the same as:

   is=x...*tol/10
...
   if abs(i1-i2)<=abs(is)

or:

   is=abs(x...)*tol/10
...
   if abs(i1-i2)<=is

[HansOlsson]

Minor note: Changing y to abs(y) ensures that the integral is non-negative, since all of the quadrature-weights are positive - (it's still a bit of a hack since abs(y) may not be sufficiently differentiable). That's normal for quadrature formulas. However, time-integration algorithms sometimes fail that (like Nyström's and Fehlberg's methods with one negative B-value each).

[AHaumer]

I've tried it: works as expected. Thanks a lot!
@HansOlsson should I add a test example?
Maybe I should add the test example in my PR #4800 after merging this one:
We can compare an analytical solution, Lobatto and different Newton-Cotes.

[AHaumer]

Well that's code hard for me to analyze, but I've tried it and it works fine!

[AHaumer]

solves #4807