[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

st: help with my simple integration command in mata

From   Nathan Danneman <>
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,

mata clear

function NDintegration(								// to perform integration via simpson's
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 luck

f = &funky()

NDintegration(f, 0, 10)


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

*   For searches and help try:

© Copyright 1996–2015 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index