Menu
▾
▴
#745 improved numerical integration in bivariat.dem
I've found a simple way to avoid the well-documented recursion limit in gnuplot's numerical integration demo, bivariat.dem. Simply change these lines:
integral_f(x) = (x>0)?int1a(x,x/ceil(x/delta)):-int1b(x,-x/ceil(-x/delta)) int1a(x,d) = (x<=d*.1) ? 0 : (int1a(x-d,d)+(f(x-d)+4*f(x-d*.5)+f(x))*d/6.) int1b(x,d) = (x>=-d*.1) ? 0 : (int1b(x+d,d)+(f(x+d)+4*f(x+d*.5)+f(x))*d/6.)
to read:
integral_f(x) = (x>0)?int1a(x,x/ceil(x/delta)):-int1b(x,-x/ceil(-x/delta)) int1a(x,d) = (x<=d*.1) ? 0 : sum [i=1:int(x/d)] ((f(i*d-d)+4*f(i*d-d*.5)+f(i*d))*d/6.) int1b(x,d) = (x>=-d*.1) ? 0 : sum [i=1:int(-x/d)] ((f(i*d+d)+4*f(i*d+d*.5)+f(i*d))*d/6.)
Could you please provide more context for this?
Is there an actual problem somewhere that you are trying to fix?
What recursion limit? What documentation?
Are you referring to this discussion from 6 years ago?
http://gnuplot.10905.n7.nabble.com/Integration-in-gnuplot-IS-possible-td11075.html
It is true that there is an impossible recursion later in "bivariat.dem" - the Ackerman function.
I think whoever wrote the demo was having a little joke, since the entire point of the Ackerman function is that it becomes impossibly expensive to compute despite its simple appearance.
Yes, that's the thread I was referring to.
The actual problem that I'm trying to fix is numerical integration of the sinc function, which exhibits visible artifacts (using the "bivariat.dem" code) when "delta" is set to .2 (the default), but triggers errors (recursion depth limit exceeded) when "delta" is decreased to smaller numbers.
I've improved the code slightly since yesterday, as it couldn't take integer arguments (only floats). This version seems to work with either integers or floats:
If I plug your new definitions (either version) into bivariat.dem with no other change, the results are ugly compared to the unmodified demo. Are you sure this is an improvement?
I have not thought deeply about this area, but it occurs to me that version 5 offers other ways to approach this problem. What if instead of setting a small delta and evaluating recursively, you set samples=10000, evaluated f(x) into an array just once, then applied Simpson's rule to the array values?
Last edit: Ethan Merritt 2017-01-19
Yes, I'm sure this is an improvement! (there's just some bugs in my code, that's all)
I found two bugs: a rounding error, and missing signs in the negative branch of the code.
Try this version:
You have an interesting idea. I made this change to the code because I was sure that there were ways to avoid the recursion limit, and this seemed like about the simplest change to do it.
OK. Now the demo output looks the same to my eye. That's good as far as it goes but doesn't show the need for better intXXX() functions. If you want to update the demo then it would be good to add your sinc(x) test case where (I gather) the original set of equations fails and your new set is works.
Even your new equations use recursion though, so I think it would be interesting to also show a non-recursive solution using an intermediate array if that is feasible.
The whole demo could use an overhaul, including changing the name to something more obvious that "bivariat". Maybe break it into smaller pieces, one of which focuses on these numerical integration plots? Maybe add a mention in the manual so that people don't have to hunt up 6 year old threads on Nabble?