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/
