[R] unexpected result in function valuation

[Ravi Varadhan]

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.