On Nov 15, 2010, at 7:18 PM, Eduardo de Oliveira Horta wrote:
Hello,
I was trying to build some functions which I would like to integrate
over an
interval using the function 'integrate' from the 'stats' package. As
an
example, please consider the function
h(u)=sin(pi*u) + sqrt(2)*sin(pi*2*u) + sqrt(3)*sin(pi*3*u) +
2*sin(pi*4*u)
Two alternative ways to 'build' this function are as in f and g below:
coeff<-sqrt(1:4)
zeta<-function(i) {force(i); function(u){
zeta<-sqrt(2)+sin(pi*i*u)
}}
f<-function(u){
f<-0
for (j in 1:4){
f<-f+coeff[j]*(zeta(j)(u))
}
f<-f
}
g<-function(u){
g<-crossprod(coeff,zeta(1:4)(u))
}
Obviously, f and g are equivalent, but in the actual code I am
writing, g is
much faster. However, I can't seem to integrate (neither to plot) g.
In
fact, these are the error messages I get:
print(f(.1))
[1] 4.443642
print(g(.1))
[,1]
[1,] 4.443642
Everything ok until here... but...
When using plot(), I get this:
plot(f) #plots, as expected.
plot(g)
Try vectorizing g()
g<-function(u){
g<-crossprod(coeff,zeta(1:4)(u))
}
gV <- Vectorize(g)
integrate(gV,0,1)
## 9.696303 with absolute error < 2.5e-13
You didn't say what you thought the analytic answer might be so I
cannot check it.
--
David.
Error in crossprod(coeff,zeta(1:4)(u)) : arguments not compatible
Besides that: Warning message:
In pi * i * u :
longer object length is not a multiple of shorter object length
When using integrate(), I get this:
integrate(f,0,1)
1.004172 with absolute error < 2.5e-13
integrate(g,0,1)
Error in crossprod(coeff,zeta(1:4)(u)) : arguments not compatible
Besides that: Warning message:
In pi * i * u :
longer object length is not a multiple of shorter object length
I have already tried some simple 'solutions', for example to set
g<-function(u){
vec1<-drop(coeff)
vec2<-drop(zeta(1:4)(u))
g<-drop(crossprod(coeff,zeta(1:4)(u)))
}
as well as using the %*% operation, but these won't solve my problem.
Any suggestions would be welcome. Thanks in advance,
Eduardo Horta
[[alternative HTML version deleted]]
David Winsemius, MD
West Hartford, CT
______________________________________________
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.
