I apologize for the noise. I didn't clean up the code enough. See below.
<<<snip>>>
>
> Craig,
>
> I haven't seen an answer to this yet, so let me jump in. You seem to have
> some stuff still leftover from MATLAB. Here is some cleaned up code that
> produces the result you expect. I don't think the value of dx was being
> correctly computed in your code. I did not change the assignment operator
> you used (=), but in R the "preferred" operator is "<-" (without the
> quotes).
>
> myquadrature <- function(f,a,b){
> npts = length(f)
> nint = npts-1
> if(npts <= 1) error('need at least two points to integrate')
> if(b <= a) error('something wrong with the interval, b should be greater
> than a') else dx=b/nint
The 2 'if' statements above should have been
if(npts <= 1) stop('need at least two points to integrate')
if(b <= a) stop('something wrong with the interval, b should be greater than
a') else dx=b/nint
> sum(f[-npts]+f[-1])/2*dx
> }
>
> #Call my quadrature
> x = seq(0,2000,10)
> h = 10*(cos(((2*pi)/2000)*(x-mean(x)))+1)
> u = (cos(((2*pi)/2000)*(x-mean(x)))+1)
> a = x[1]
> b = x[length(x)]
> plot(x,-h)
> a = x[1];
> b = x[length(x)]
>
> #call your quadrature function. Hint, the answer should be 30000.
> f = u*h
> result = myquadrature(f,a,b)
> result
>
> Hope this is helpful,
>
> Dan
>
Daniel Nordlund
Bothell, WA USA
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
