# Numerical Integration through Plotting

Inspired by this on Spiked Math:

(source: spikedmath.com)

In this code golf, you are going to write a numerical integrator. Instead of weighting, you first plot on a graphic canvas pixel-by-pixel or plot on a 2D array, then calculate the integral by counting the colored pixels.

You are implementing a function with the following arguments:

• f: the function, in lambda or function pointer
• xmin, xmax: the interval, also the familiar window settings on a typical graphic calculator
• ymin, ymax: window settings for y-axis, note: it is the caller's responsibility to choose ymax, because the function can go over ymax, and the parts went over will not be counted.
• xres, yres: resolution of the plot: size of the canvas or the dimensions of the 2D array

Here is an explanation:

To integrate sin(x) from 1.0 to 8.0 with y = -1.0 to 1.0, and resolution 30 * 15, first plot

 ++++                       ++
++++++                     +++
+++++++                   ++++
++++++++                 +++++
++++++++                 +++++
+++++++++               ++++++
+++++++++              +++++++
-------------
-----------
-----------
---------
--------
-------
-----


Count the total number of +s: 83 total number of -s: 64 and N = 83 - 64 = 19

(the formula for N is actually notice xres-1, x never equal to xres, don't make this off-by-one error )

and calculate the integral:

N / canvas size * window area
= (N / (xres * yres) ) * ((ymax - ymin) * (xmax - xmin))
= 19 / (30 * 15) * ((1 - (-1)) * (8 - 1))
~= 0.5911


Shortest code wins.

Test cases:

| f(x)       |  xmin |  xmax |  ymin | ymax  | xres | yres | S f(x) dx |
| 3          |  0.0  |  10.0 |   0.0 |   2.0 |   25 |   25 |  20.0     |
| 3          |  0.0  |  10.0 |   0.0 |   5.0 |    5 |    5 |  30.0     |
| 3          |  0.0  |  10.0 |   0.0 |   4.0 |    5 |    5 |  32.0     |
| 2x + 2     | -5.0  |   5.0 | -10.0 |  10.0 |   50 |   50 |  18.72    |
| 2x + 2     | -5.0  |   5.0 | -20.0 |  20.0 |  500 |  500 |  20.1648  |
| x^2        |  0.0  |  10.0 |   0.0 | 120.0 |   10 |   10 | 348.0     |
| x^2        |  0.0  |  10.0 |   0.0 | 120.0 |  100 |  100 | 335.88    |
| x^2        |  0.0  |  10.0 |   0.0 | 120.0 |  500 |  500 | 333.6768  |
| x^2 sin(x) | -3.14 |  3.14 |  -5.0 |   5.0 |   50 |   50 |   0.60288 |
| x^2 sin(x) | -3.14 |  3.14 |  -5.0 |   5.0 | 1500 | 1500 |   0.60288 |
| sin(x^2)   |  0.0  | 10.0  |  -1.0 |   1.0 |  250 |   50 |   0.7808  |
| sin(x^2)   |  0.0  | 10.0  |  -1.0 |   1.0 | 1000 | 1000 |   0.59628 |

• This is not well specified. Will the function be supplied as a function argument? Or is the code supposed to parse some expression syntax? Is the aim a function, or a program that takes text input and output? Does the code need to actually plot into a 2D array, or just calculate the number of pixels that would be turned on? If the output of the plot is optional, then no one will do it. Apr 17 '11 at 3:38
• So quadrature by the midpoint rule? Apr 17 '11 at 6:53
• no, integrate by counting pixels on the plot, something like (blue pixels - red pixels) / (total pixels) * (xmax - xmin) * (ymax - ymin) Apr 17 '11 at 18:07
• Ok, so it's left point rather than midpoint and quantised. What correctness constraints are there? I'm sure I could construct a test case which distinguishes different ways of looping over the x-values. Apr 17 '11 at 22:37
• It's Monte Carlo without randomness :D Apr 18 '11 at 1:45

function I(f,x,X,y,Y,r,R)I=io.write N=0 for k=R-1,0,-1 do for l=0,r-1 do v=f(x+l*(X-x)/r) K=y+k*(Y-y)/R c=(K<0 and K>v)and'-'or(K>0 and K<v and'+')or' 'N=c=="+"and N+1 or(c=="-"and N-1 or N)I(c)end I'\n'end return N end

assert(I(math.sin,1.0,8.0,-1,1,30,15)==19)