Math.Nonlinear.quadratureNewtonCotes

[AHaumer]

When investigating the relation between spatial flux density distribution of a rotating rotor (with permanent magnets) and the flux linkage of a coil in the stator core (to calculate the induce voltage), I have to calculate an integral over the angle spanned by the coil. Since quadratureLobatto has some problems with a moving integration interval, I implemented two less sophisticated quadrature functions: quadratureTrapezoidal (based on the trapezoidal rule) and quadratureSimpson (based on Simpson's rule).
Both are tested in the example MovingIntegral (using the utility function funB). That's also a nice example for the usage of quadrature functions. Defined comparison signals, no backwards compatibility issues.

[AHaumer]

I managed to implement the first 4 closed Newton-Cotes formulas in 1 function, just distinguished by 1 parameter.
This is the desired formulation, allowing the user to choose the degree of the interpolation polynominal and the number of intervals.

[AHaumer]

pinging potential reviewers @christiankral @casella @HansOlsson @henrikt-ma @maltelenz @MartinOtter

[AHaumer]

@HansOlsson to explain the motivation to implement rather simple quadrature formulas:
I tried to use Modelica.Math.Nonlinear.quadratureLobatto, but all 3 tools I've used for testing get stuck (see enclosed model).
The Newton-Cotes-formulas in this PR work fine on this "moving integral" with the same 3 tools.

MovingIntegralSinusoidal.zip

[HansOlsson]

Looking at the code and at: https://www.researchgate.net/publication/226706221_Adaptive_Quadrature-Revisited (which should replace the ftp-link) my only reaction is that the code is stupid and stupidly obscure.

Without looking further: Do you see and understand the stopping criteria for the recursion, and how it relates to tolerance?

The correct answer is:
if(is+(i1-i2)==is)
where is is computed as original integral times tolerance/eps!

There are a number of issues with that:

• It is really really fragile
  • Real comparisons are inherently fragile,
  • Here it adds the logic that (i1-i2) should be evaluated first (maybe happens in Matlab, not guaranteed in Modelica) and the addition not computed with extra bits (may happen).
• As a pure relative tolerance it fails badly if the integral is basically zero. That happens at time=0.25.
• Probably more

[AHaumer]

@HansOlsson no I don't understand the code, and I don't trust it.
That's why I implemented an alternative much simpler solution.
I do not know who implemented quadratureLobatto. Main authors of Math: Martin Otter and Marcus Baur.
@MartinOtter could you please comment on the issues with quadratureLobatto?

[HansOlsson]

@HansOlsson no I don't understand the code, and I don't trust it. That's why I implemented an alternative much simpler solution. I do not know who implemented quadratureLobatto. Main authors of Math: Martin Otter and Marcus Baur. @MartinOtter could you please comment on the issues with quadratureLobatto?

As far as I understand the code implements the algorithm from the paper, it's just that the paper is more trying to be clever than correct and robust; so the problem isn't so much the Modelica-implementation as the underlying paper (although there are at least two issues with the Modelica implementation - related to the convergence criteria).

Looking more at this and the paper contents:

• The authors are proud of their is+(i1-i2)==is trick!
• They add a guard against is being exactly zero, without considering what happens if is is very close to zero.
• The Modelica code uses 10*Modelica.Constans.eps instead of Modelica.Constans.eps. One might think that it loosens up tolerances, but as far as I can tell the consequence is that effectively the user given tolerances are divided by a factor of 10 (even if they are 1e-6 or 1e-7).

I have tried to patch it as follows:
Use abs(i1-i2)<is instead of is+(i1-i2)==is, and skip division by eps when computing is.
Replace tolerance formula by is := max(abs(is),isabs/10)*tol/10; where isabs is computed using abs of is-formula; to guard against almost zero values.