6
\$\begingroup\$

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 formula 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 |
\$\endgroup\$
  • 1
    \$\begingroup\$ 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. \$\endgroup\$ – MtnViewMark Apr 17 '11 at 3:38
  • \$\begingroup\$ So quadrature by the midpoint rule? \$\endgroup\$ – Peter Taylor Apr 17 '11 at 6:53
  • \$\begingroup\$ no, integrate by counting pixels on the plot, something like (blue pixels - red pixels) / (total pixels) * (xmax - xmin) * (ymax - ymin) \$\endgroup\$ – Ming-Tang Apr 17 '11 at 18:07
  • \$\begingroup\$ 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. \$\endgroup\$ – Peter Taylor Apr 17 '11 at 22:37
  • \$\begingroup\$ It's Monte Carlo without randomness :D \$\endgroup\$ – Eelvex Apr 18 '11 at 1:45
2
\$\begingroup\$

Lua, 219 chars

Done it all in 1 step, no need for an intermediate table. Plots the function too ;)

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

Test code:

assert(I(math.sin,1.0,8.0,-1,1,30,15)==19)
| improve this answer | |
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.