Dear all,
I have a very small script to plot a function. Here it is:
##########################################
sinca <- function(N,th)
{
return(sin((N+0.5)*th)/sin(0.5*th))
}
plot_sinca <- function(N)
{
x <- seq(-5*pi,5*pi,by=pi/100)
y <- rep(0,length=length(x))
for (i in 1:length(x))y[i] <- sinca(N,x[i])
plot(x,y,type="l",ylim=c(0,2*N+4))
return(c(x,y))
}
##########################################
When I load the script and run the function like this:
###########################################> data <- plot_sinca(4)
> y <- data[1002:2002]
###########################################
I notice a spike on the plot which should not be there.
In fact I have checked and:
###########################################> y[701]
[1] 10.07404> sinca(4,2*pi)
[1] 9
###########################################
The second result is the correct one. Why, then do
I get the y[701]=10.07404? This function is not supposed
to be higher than 9...
Any help is greatly appreciated.
Regards,
J
Dr James Foadi
Membrane Protein Laboratory
Diamond Light Source Ltd
Chilton, Didcot
Oxfordshire OX11 0DE
---
[[alternative HTML version deleted]]
What value should your formula give when x is a multiple of 2*pi? You seem to believe 9 is correct but in fact NaN is. Element 701 of x is approximately but not exactly 2*pi: on my system it is about 7*.Machine$double.eps different. You cannot expect sin(N*pi) to be exactly zero for N != 0. On Thu, 5 Jul 2007, James Foadi wrote:> Dear all, > I have a very small script to plot a function. Here it is: > > ########################################## > sinca <- function(N,th) > > { > > return(sin((N+0.5)*th)/sin(0.5*th)) > > } > > plot_sinca <- function(N) > > { > > x <- seq(-5*pi,5*pi,by=pi/100) > > y <- rep(0,length=length(x)) > > for (i in 1:length(x))y[i] <- sinca(N,x[i]) > > plot(x,y,type="l",ylim=c(0,2*N+4)) > > return(c(x,y)) > > } > > ########################################## > > When I load the script and run the function like this: > > ########################################### >> data <- plot_sinca(4) >> y <- data[1002:2002] > ########################################### > > I notice a spike on the plot which should not be there. > In fact I have checked and: > ########################################### >> y[701] > [1] 10.07404 >> sinca(4,2*pi) > [1] 9 > ########################################### > > The second result is the correct one. Why, then do > I get the y[701]=10.07404? This function is not supposed > to be higher than 9... > > Any help is greatly appreciated. > > Regards, > > J > > Dr James Foadi > Membrane Protein Laboratory > Diamond Light Source Ltd > Chilton, Didcot > Oxfordshire OX11 0DE > --- > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at stat.math.ethz.ch 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. >-- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595
Very simple; it is your function. You need to step through and see that you are evaluating close to zero:> x[701][1] 6.283185> sin((4.5*x[701]))[1] -1.666142e-14> sin(.5*x[701])[1] -1.653896e-15> sin((4.5*x[701]))/sin(.5*x[701])[1] 10.07404>With numbers that small you might be losing significance. See the FAQ on floating point numbers. On 7/5/07, James Foadi <james.foadi@diamond.ac.uk> wrote:> > Dear all, > I have a very small script to plot a function. Here it is: > > ########################################## > sinca <- function(N,th) > > { > > return(sin((N+0.5)*th)/sin(0.5*th)) > > } > > plot_sinca <- function(N) > > { > > x <- seq(-5*pi,5*pi,by=pi/100) > > y <- rep(0,length=length(x)) > > for (i in 1:length(x))y[i] <- sinca(N,x[i]) > > plot(x,y,type="l",ylim=c(0,2*N+4)) > > return(c(x,y)) > > } > > ########################################## > > When I load the script and run the function like this: > > ########################################### > > data <- plot_sinca(4) > > y <- data[1002:2002] > ########################################### > > I notice a spike on the plot which should not be there. > In fact I have checked and: > ########################################### > > y[701] > [1] 10.07404 > > sinca(4,2*pi) > [1] 9 > ########################################### > > The second result is the correct one. Why, then do > I get the y[701]=10.07404? This function is not supposed > to be higher than 9... > > Any help is greatly appreciated. > > Regards, > > J > > Dr James Foadi > Membrane Protein Laboratory > Diamond Light Source Ltd > Chilton, Didcot > Oxfordshire OX11 0DE > --- > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help@stat.math.ethz.ch 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. >-- Jim Holtman Cincinnati, OH +1 513 646 9390 What is the problem you are trying to solve? [[alternative HTML version deleted]]
The problem is that you are dividing two numbers that are both very small.
Any small imprecision in the denominator creates a big error.
Here you can re-write your function using a trig identity to get accurate
results:
sinca <- function(N,th) {
#return(sin((N+0.5)* th)/sin(0.5*th))
return( (sin(N*th)/tan(th/2)) + cos(N*th))
}
This function works well, as you had expected.
Ravi.
----------------------------------------------------------------------------
-------
Ravi Varadhan, Ph.D.
Assistant Professor, The Center on Aging and Health
Division of Geriatric Medicine and Gerontology
Johns Hopkins University
Ph: (410) 502-2619
Fax: (410) 614-9625
Email: rvaradhan at jhmi.edu
Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html
----------------------------------------------------------------------------
--------
-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of James Foadi
Sent: Thursday, July 05, 2007 1:39 PM
To: r-help at stat.math.ethz.ch
Subject: [R] unexpected result in function valuation
Dear all,
I have a very small script to plot a function. Here it is:
##########################################
sinca <- function(N,th)
{
return(sin((N+0.5)*th)/sin(0.5*th))
}
plot_sinca <- function(N)
{
x <- seq(-5*pi,5*pi,by=pi/100)
y <- rep(0,length=length(x))
for (i in 1:length(x))y[i] <- sinca(N,x[i])
plot(x,y,type="l",ylim=c(0,2*N+4))
return(c(x,y))
}
##########################################
When I load the script and run the function like this:
###########################################> data <- plot_sinca(4)
> y <- data[1002:2002]
###########################################
I notice a spike on the plot which should not be there.
In fact I have checked and:
###########################################> y[701]
[1] 10.07404> sinca(4,2*pi)
[1] 9
###########################################
The second result is the correct one. Why, then do
I get the y[701]=10.07404? This function is not supposed
to be higher than 9...
Any help is greatly appreciated.
Regards,
J
Dr James Foadi
Membrane Protein Laboratory
Diamond Light Source Ltd
Chilton, Didcot
Oxfordshire OX11 0DE
---
[[alternative HTML version deleted]]
______________________________________________
R-help at stat.math.ethz.ch 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.