[R] MATLAB vrs. R

[Daniel Nordlund]

> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
> On Behalf Of Craig O'Connell
> Sent: Monday, October 11, 2010 8:10 AM
> To: alain.guillet at uclouvain.be
> Cc: r-help at r-project.org; pdalgd at gmail.com
> Subject: Re: [R] MATLAB vrs. R
> 
> 
> alain,
> 
>     Perhaps i'm still entering the code wrong.  I tried using your
> result=myquadrature(f,0,2000)
> print(result)
> 
> Instead of my:
> val = myquadrature(f,a,b)
> result=myquadrature(val,0,2000)
> print(result)
> 
> ...and I am still getting an inf inf inf inf inf...
> 
> Did you change any of the previous syntax in addition to changing the
> result statement?
> 
> Thank you so much and I think my brain is fried!  Happy Holiday.
> 
> Craig
> 

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