# st: help with my simple integration command in mata

 From Nathan Danneman To statalist@hsphsun2.harvard.edu Subject st: help with my simple integration command in mata Date Mon, 15 Jun 2009 15:07:15 -0400

```Hi all,

I have written the following bit of code that uses Simpson's method to
approximate definite integrals.  It doesn't seem to work: it responds
that there is a conformability error, which I cannot seem to locate.

I would welcome any general tips on debugging programs in mata, and
especially welcome any specific tips on refining this code.

Thanks so much,
Nathan

mata
mata clear

function NDintegration(								// to perform integration via simpson's
approximation
pointer(real scalar function) scalar f,						//points to the function
to be integrated
real scalar lower,									//lower bound of integration
real scalar higher)									//upper bound
{

if (lower > higher) return("limits of integration reversed")
//make sure upper > lower

width = (higher-lower) / 14								//calculates width of each "slice"

Y = (.)											//creates an empty matrix

for (i=0; i<=14; i++) {
x = lower + (i*width)
yi = (*f)(x)								//is this part correct?  I'm not sure I've used
the pointer correctly
d = yi
Y = (Y, d)
}
st_matrix("Y", Y)

return( (width/3) * (1*st_matrix("Y")[1,2] + 2*st_matrix("Y")[1,3] +
4*st_matrix("Y")[1,4] + 2*st_matrix("Y")[1,5] + 4*st_matrix("Y")[1,6]
+ 2*st_matrix("Y")[1,7] + 4*st_matrix("Y")[1,8] +
2*st_matrix("Y")[1,9] + 4*st_matrix("Y")[1,10] +
2*st_matrix("Y")[1,11] + 4*st_matrix("Y")[1,12] +
2*st_matrix("Y")[1,13] + 4*st_matrix("Y")[1,14] +
1*st_matrix("Y")[1,15]) )

//this last bit simply sums up the y values and multiplies them by the
appropriate fudge factors

}

function funky(real scalar x) {                       //I tried to use
this simple function to test the integration routine...no luck
10*x
}

f = &funky()

NDintegration(f, 0, 10)

end

--
"Imagination is more important than knowledge..."
-- Albert Einstein

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/
```