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);
Thank you and any potential suggestions would be greatly appreciated.
Dr. Argese.
[[alternative HTML version deleted]]
______________________________________________
[email protected] 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.
