Re: [R] MATLAB vrs. R
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 3.
f=u.*h;
val = myquadrature(f,a,b);
fprintf('the solution is %f\n',val);
This is great, I got the expected answer of 3.
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.
Re: [R] MATLAB vrs. R
Thank you Peter. That is very much helpful. If you don't mind, I continued
running the code to attempt to get my answer and I continue to get inf inf
inf... (printed around 100 times).
Any assistance with this issue. Here is my code (including your corrections):
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)
int=0
int - 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 = 1.*(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 3.
f=u*h;
val = myquadrature(f,a,b); ? ___This is where issue arises.
result=myquadrature(val,0,2000) ?
print(result) ?
Thanks again,
Phil
[[alternative HTML version deleted]]
__
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.
Re: [R] MATLAB vrs. R
Hi,
The first argument of myquadrature in result shouldn't be val but f I
guess. At least it works for me
result=myquadrature(f,0,2000)
print(result)
[1] 3
Regards,
Alain
On 11-Oct-10 09:37, Craig O'Connell wrote:
Thank you Peter. That is very much helpful. If you don't mind, I continued
running the code to attempt to get my answer and I continue to get inf inf
inf... (printed around 100 times).
Any assistance with this issue. Here is my code (including your corrections):
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)
int=0
int- 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 = 1.*(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 3.
f=u*h;
val = myquadrature(f,a,b); ? ___This is where issue arises.
result=myquadrature(val,0,2000) ?
print(result) ?
Thanks again,
Phil
[[alternative HTML version deleted]]
__
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.
--
Alain Guillet
Statistician and Computer Scientist
SMCS - IMMAQ - Université catholique de Louvain
Bureau c.316
Voie du Roman Pays, 20
B-1348 Louvain-la-Neuve
Belgium
tel: +32 10 47 30 50
__
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.
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
Date: Mon, 11 Oct 2010 09:59:17 +0200
From: alain.guil...@uclouvain.be
To: craigpoconn...@hotmail.com
CC: pda...@gmail.com; r-help@r-project.org
Subject: Re: [R] MATLAB vrs. R
Hi,
The first argument of myquadrature in result shouldn't be val but f I
guess. At least it works for me
result=myquadrature(f,0,2000)
print(result)
[1] 3
Regards,
Alain
On 11-Oct-10 09:37, Craig O'Connell wrote:
Thank you Peter. That is very much helpful. If you don't mind, I continued
running the code to attempt to get my answer and I continue to get inf inf
inf... (printed around 100 times).
Any assistance with this issue. Here is my code (including your
corrections):
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)
int=0
int- 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 = 1.*(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 3.
f=u*h;
val = myquadrature(f,a,b); ? ___This is where issue arises.
result=myquadrature(val,0,2000) ?
print(result) ?
Thanks again,
Phil
[[alternative HTML version deleted]]
__
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.
--
Alain Guillet
Statistician and Computer Scientist
SMCS - IMMAQ - Université catholique de Louvain
Bureau c.316
Voie du Roman Pays, 20
B-1348 Louvain-la-Neuve
Belgium
tel: +32 10 47 30 50
[[alternative HTML version deleted]]
__
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.
Re: [R] MATLAB vrs. R
-Original Message-
From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org]
On Behalf Of Craig O'Connell
Sent: Monday, October 11, 2010 8:10 AM
To: alain.guil...@uclouvain.be
Cc: r-help@r-project.org; pda...@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 3.
f = u*h
result = myquadrature(f,a,b)
result
Hope this is helpful,
Dan
Daniel Nordlund
Bothell, WA USA
__
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.
Re: [R] MATLAB vrs. R
I apologize for the noise. I didn't clean up the code enough. See below.
snip
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
The 2 'if' statements above should have been
if(npts = 1) stop('need at least two points to integrate')
if(b = a) stop('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 3.
f = u*h
result = myquadrature(f,a,b)
result
Hope this is helpful,
Dan
Daniel Nordlund
Bothell, WA USA
__
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.
Re: [R] MATLAB vrs. R
Daniel,
That's it! Thanks. Your help is very much appreciated. I'm hoping to
nail down the code conversion from MATLAB to R, but it seems to be a bit more
difficult that I had anticipated.
Craig
From: djnordl...@frontier.com
To: djnordl...@frontier.com; r-help@r-project.org
Date: Mon, 11 Oct 2010 16:12:48 -0700
Subject: Re: [R] MATLAB vrs. R
I apologize for the noise. I didn't clean up the code enough. See below.
snip
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
The 2 'if' statements above should have been
if(npts = 1) stop('need at least two points to integrate')
if(b = a) stop('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 3.
f = u*h
result = myquadrature(f,a,b)
result
Hope this is helpful,
Dan
Daniel Nordlund
Bothell, WA USA
__
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.
[[alternative HTML version deleted]]
__
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.
[R] MATLAB vrs. R
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 3.
f=u.*h;
val = myquadrature(f,a,b);
fprintf('the solution is %f\n',val);
This is great, I got the expected answer of 3.
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]]
__
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.
