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....@cbs.dk Priv: pda...@gmail.com
______________________________________________
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.
