Problems with Romberg and Simpson integration

[marcobonici]

Dear developer, hi to everybody! I am new to Julia, I am physicist and I am trying to understand if it makes sense to use julia for my next scientific projects. I was trying to perform some numerical integration and I found your package, so I decided to make some simple tests.
k = 2^4
x = LinRange(0, π, k+1)
y = sin.(x)
println(integrate(x, y, Trapezoidal()))
println(integrate(x, y, SimpsonEven()))
println(integrate(x, y, RombergEven()))

obtaining
1.9935703437723393
2.222549924725271
2.000005549979671
In order to be sure I am doing things correctly, I decided to cross-check the results performing the same calculation in Python
k = 2**4
x = np.linspace(0, np.pi, k+1)
y = np.sin(x)
print(np.trapz(y,x))
print(integrate.simps(y,x))
print(integrate.romb(y=y, dx=x[1]-x[0], show=False))
obtaining
1.9935703437723393
2.0000165910479355
1.9999999945872902
As you can see, the trapz integrations looks fine but there are huge discrepancies in the Simpson and Romberg results. I looked into your issues and I found this Romber.jl, so I decided to use it finding
println(romberg(x, y))
(1.9999999945872906, 5.555392379896773e-6)
This result is almost the same of the scipy simpson method!
Furthermore, I found that if try
k = 2^5
x = LinRange(0, π, k+1)
y = sin.(x)
println(integrate(x, y, RombergEven()))
I obtain
1.9999999945872906
Which is exactly the result of Romberg.jl and scipy! it looks that the result of RombergEven with 2^(k+1)+1 points is the same of scipy and romberg with 2^k+1 points.
Am I using your package in the wrong way or there is something else?
Thank you for your time!

[dextorious]

Thank you for the report. The Simpson issue you're having has been resolved for a while, but the version that has the fixes hasn't yet been upstreamed into JuliaRegistry (that's being dealt with in an adjacent issue).

Running your example on the master branch yields:

julia> k = 2^4; x = LinRange(0, π, k+1); y = sin.(x);

julia> integrate(x, y, Trapezoidal())
1.9935703437723395

julia> integrate(x, y, SimpsonEven())
2.0000183531747346

julia> integrate(x, y, RombergEven())
┌ Warning: RombergEven :: final step reached, but accuracy not: 0.001434818155835238 > 1.0e-12
└ @ NumericalIntegration C:\Users\dxt\.julia\dev\NumericalIntegration\src\NumericalIntegration.jl:182
2.000005549979671

The Romberg integration still has an accuracy issue, but it should reliably output warnings on the edge cases where it cannot reach the expected accuracy due to the algorithm used. This is a compromise between performance and accuracy, since historically this package was basically written to get the fastest performance possible. In the past few days Romberg.jl also got a new implementation that might be able to compete on performance as well (also described in a fellow issue here), so I might just add that as a fallback dependency to cover these edge cases as well.

For the moment your options are either running the master branch (]add NumericalIntegration#master) or waiting until the JuliaRegistry / version tagging issue is resolved, which should hopefully happen soon, as well as using Romberg.jl if guaranteed accuracy is a priority. Either way, I'm happy to welcome a fellow physicist to the community. 👍

[marcobonici]

Hi @dextorious ,
Thank you so much for your answer! As I said, I just started studying Julia, so I am just wandering around:)
Just a few more questions:

1. The method SimpsonEven you implemented is the alternative simpson method (as you have written inside the code), which confused me at the beginning: when I tried the function I was expecting to find the same result of scipy (which implements "classic" simpson rule). This is a matter of choice (and yours is perfectly legit) but maybe naming the function AlternativeSimpsonEven would be more transparent to users? Probably Python user will make a comparison with scipy's method as I have done and the result will be different since the algorithm implemented is actually different, not due to an issue.
2. This more a curiosity since I am new to Julia. I see that at the beginning of the package you defined an abstract type, IntegrationMethod, and some subtypes. The different integration methods are invoked using the third argument of the integrate function. This is an application of multiple dispatch and the most "Julia way" to do this kind of thing, right? Sorry to bother you, but I'd like to be sure I have understood correctly this item, since I am going to use it soon:)
   Best regards,
   Marco

[dextorious]

1. To be honest, I'm not sufficiently familiar with the literature to clarify the nomenclature. I gave the wikipedia reference for clarity, but I've seen both versions referred to as just "Simpson's rule" before in textbooks (after all, the rule is the same - the difference is in which segments you apply it to). I think both yield roughly comparable accuracy, but I certainly don't give any guarantees of identical results with any other implementation, SciPy's or anyone else's. If there's any legitimate advantage to SciPy's version, I'm happy to accept a pull request that adds it as well, but generally this package aims for simple and fast implementations over covering every possible method or edge case. There's already a number of packages that aim for more sophisticated methods and higher accuracy at the cost of lower performance or more complicated APIs, so there's no need to do everything here.

2. Julia doesn't really enforce a single paradigm in any way, but it is Julia-specific in the sense that it relies on multiple dispatch, meaning there aren't many languages out there that support that kind of approach. The advantage here is performance - some edge cases aside, the dispatch is optimized away and there's no runtime cost like there would be with a runtime branch (like an if statement checking which method to use). Whether that makes the API better/clearer or not is somewhat subjective, though many in the Julia community seem inclined to prefer it this way, myself included.