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.

[HansOlsson]

I agree that the original adaptive quadratureLobatto wasn't working - and needs to be updated; see separate PR.

However, a non-adaptive Newton-Cotes as this one seems like a step backwards, so:

• It may be a necessary step backwards
• Making it adaptive seem possible, but might run into the same issues and clearly take a lot of time to get right
• A compromise might be to add an optional output for an error estimate
• Clearly state that it is non-adaptive and highlight the benefits of a non-adaptive algorithm, e.g., that the result is differentiable

The last option may be the best.

[HansOlsson]

Looks ok, but might better document that it is non-adaptive which gives differentiable results.

[AHaumer]

@HansOlsson differentiable with respect to what? fun? a and b?
At least in my example it doesn't work to use v=-der(psi);

Error term: Good idea but rather hard to compute, you need a higher derivative of the integrand (taken from Wikipedia):

The variable ξ should be varied over the whole interval [a, b] and pick the highest value of the derivative.

So the only thing that remains is to state that the formulas are non-adaptive but "composite" fomrulations,
i.e. dividing the interval [a, b] into sub-intervals and apply the Newton-Cotes formula to each sub-interval.