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
with your program though. Your loop calculates 15 values (0 through
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
readable
}
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 <[email protected]> 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/
