RE: st: help with my simple integration command in mata

 From Glenn Goldsmith To statalist@hsphsun2.harvard.edu Subject RE: st: help with my simple integration command in mata Date Tue, 16 Jun 2009 00:36:52 +0100

Hi Nathan,

1. The error is occurring because your test function -funky()- needs
to -return()- 10*x rather than just displaying it as it does now.
Replacing the -10*x- line with -return(10*x)- should do the trick.

2. It looks like there might be something else slightly odd going on
14) but your return statement only uses 14 of these (2 through 15). I
assume it should use all 15, with the multipliers
(1,4,2,4,2,...,2,4,1) rather than the (1,2,4,2,...,2,4,1) you
currently have?

3. In terms of general debugging, mata can be a little difficult; but
inserting -mata set matalnum on- after your -mata clear- line will at
least give you line numbers for the errors. (You'll want to turn it
off after testing though, because when it's on it stops the compiler
from optimizing the code.)

4. There are a couple of other ways you could clean up your code to
achieve what you want slightly more efficiently. A modified version is
below:

**********************************************************************************
mata
mata clear
mata set matalnum off

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

real scalar width
//I like declarations, your mileage may vary
real colvector Y
real rowvector X

if (lower > higher) _error("limits of integration reversed")		//make
sure upper > lower NB: better to use _error() than return()

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

Y = J(15,1,.)										//creates an empty 15 x 1 column vector
for (i=1; i<=15; i++) {
x = lower + ((i-1)*width)
Y[i,1]  = (*f)(x)								        //places function evaluations
directly into Y, avoiding all the copying you were doing previously
}

X = (1,4,2,4,2,4,2,4,2,4,2,4,2,4,1)
//amended 1 x 15 row vector of multipliers
return(X*Y*width/3)
//matrix multiplication makes this quicker and more

}

scalar funky(real scalar x) {
//declares the output to be a scalar. otherwise you
could run into problems in the main program
return(10*x)
}

f = &funky()

NDintegration(f, 0, 10)

end
***********************************************************

HTH,

Glenn.

Nathan Danneman <ndanneman@gmail.com> wrote:

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
*
*   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/