[R] Trapezoid Rule

craig at sharkdefense.com craig at sharkdefense.com
Fri Oct 8 20:26:10 CEST 2010


Dear R Users,

    I've never used R before and my professor has asked us to do some
pretty intense programming (or it's intense to me at least).  Here is
the question:  Modify the function myquadrature inside the script so
that it returns the quadrature of descrete data using the trapezoidal
rule.  Modify the call to the function at the bottom of the script so
that is uses your modifies routine to compute the total volume flux
through a transect using descrete observations of bathymetry and
vertically-averaged velocity.

The professor supplied us with some sample code and i'm pretty sure it
outlines the whole problem.  My issue is, I have no idea how to write the
sum code for the trapezoid rule to estimate the area under the curve:
Here is the code provided:

myquadrature = function(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 = length(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;



# trapezoidal rule
# can code in line, hint:  sum of f is sum(f)
# f[1] is the first entry in f


}

# velocity profile across a channel
# remember to use ? for help, e.g. ?seq
x = seq(0,2000,10)
# 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[length(x)]
# this works because length(x) returns the number of values in the series
# the function cos is cosine and mean gives the mean value
# pi is 3.1415, or pi
h = 10.*(cos(((2*pi)/2000)*(x-mean(x)))+1)    #depth
u = 1.*(cos(((2*pi)/2000)*(x-mean(x)))+1)     #vertically-averaged
cross-transect velocity

# set begin and end points for the integration
a = x[1]
b = x[length(x)]
# call your quadrature function.  Hint, the answer should be 30000.
print(myquadrature())

Any assistance with this problem would be greatly appreciated.

Thanks,

Casey



More information about the R-help mailing list