[R] MATLAB vrs. R

[Peter Dalgaard]

On 10/11/2010 07:46 AM, Craig O'Connell wrote:
> 
> I need to find the area under a trapezoid for a research-related project.  I was able to find the area under the trapezoid in MATLAB using the code:
> function [int] = myquadrature(f,a,b) 
> % user-defined quadrature function
> % integrate data f from x=a to x=b assuming f is equally spaced over the interval 
> % use type 
> % determine number of data points
> npts = prod(size(f));
> nint = npts -1; %number of intervals
> if(npts <=1)
>   error('need at least two points to integrate')
> end;
> % set the grid spacing
> if(b <=a)
>   error('something wrong with the interval, b should be greater than a')
> else
>   dx = b/real(nint);
> end;
> npts = prod(size(f));
> 
> % trapezoidal rule
> % can code in line, hint:  sum of f is sum(f) 
> % last value of f is f(end), first value is f(1)
> % code below
> int=0;
> for i=1:(nint)
>     %F(i)=dx*((f(i)+f(i+1))/2);
>     int=int+dx*((f(i)+f(i+1))/2);
> end
> %int=sum(F);
> 
> Then to call "myquadrature" I did:
> % example function call test the user-defined myquadrature function
> % setup some data
> 
> % velocity profile across a channel
> % remember to use ? for help, e.g. ?seq 
> x = 0:10:2000; 
> % you can access one element of a list of values using brackets
> % x(1) is the first x value, x(2), the 2nd, etc.
> % if you want the last value, a trick is x(end)  
> % the function cos is cosin and mean gives the mean value
> % pi is 3.1415, or pi
> % another hint, if you want to multiple two series of numbers together
> % for example c = a*b where c(1) = a(1)*b(1), c(2) = a(2)*b(2), etc.
> % you must tell Matlab you want "element by element" multiplication
> % e.g.:    c = a.*b    
> % note the "."
> %
> h = 10.*(cos(((2*pi)/2000)*(x-mean(x)))+1); %bathymetry
> u = 1.*(cos(((2*pi)/2000)*(x-mean(x)))+1);  %vertically-averaged cross-transect velocity
> plot(x,-h)
> % set begin and end points for the integration 
> a = x(1);
> b = x(end);
> % call your quadrature function.  Hint, the answer should be 30000.
> f=u.*h;
> val =  myquadrature(f,a,b);
> fprintf('the solution is %f\n',val);
> 
> This is great, I got the expected answer of 30000.
>  
> NOW THE ISSUE IS, I HAVE NO IDEA HOW THIS CODE TRANSLATES TO R.  Here is what I attempted to do, and with error messages, I can tell i'm doing something wrong:
>   myquadrature<-function(f,a,b){
> npts=length(f)
> nint=npts-1
> if(npts<=1)
> error('need at least two points to integrate')
> end;
> if(b<=a)
> error('something wrong with the interval, b should be greater than a')
> else
> dx=b/real(nint)
> end;
> npts=length(f)
> _________(below this line, I cannot code)
> int=0
> for(i in 1:(npts-1))
> sum(f)=((b-a)/(2*length(f)))*(0.5*f[i]+f[i+1]+f[length(f)])}
>     %F(i)=dx*((f(i)+f(i+1))/2);
>     int=int+dx*((f(i)+f(i+1))/2);
> end
> %int=sum(F);
>  

For a literal translation, just pay a little more attention to detail:

for(i in 1:(npts-1))
	int <- int+dx*(f[1]+f[i+1])/2

However, a more R-ish way is to drop the loop and vectorize:

int <- sum(f[-npts]+f[-1])/2*dx

(or int <- sum(f) - (f[1]+f[npts])/2, by a well-known rewrite of the
trapezoidal rule).


> Thank you and any potential suggestions would be greatly appreciated.
>  
> Dr. Argese.   		 	   		  
> 	[[alternative HTML version deleted]]



-- 
Peter Dalgaard
Center for Statistics, Copenhagen Business School
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com