Thanks a lot William. I'm sorry about the syntax problem. I was working at the same time with R code and a document, so I copied from the wrong document.
I saw that I had a problem with the conditions of the integrand. A ">" instead of a ">=" on the second one. Again, thanks. One last question: Is there a way to use "approx" as the integrand? Best regards. Julio --- El vie 18-dic-09, William Dunlap <wdun...@tibco.com> escribió: > De: William Dunlap <wdun...@tibco.com> > Asunto: RE: [R] Numerical Integration > A: "Julio Rojas" <jcredbe...@ymail.com> > Fecha: viernes, 18 diciembre, 2009, 4:03 pm > The immediate problem is that > integrand(0.5) is NULL > because none of the if clauses is TRUE if x==fx[2]. > > If you restructured the code such that it had no return > statements and instead assigned things to a variable > and returned that variable at the end it would be more > apparent what happened: you would get a 'variable not > found' error when returning. > i <- function(x) { > if (x<fx[1]) retval <- 0 > else if (x>fx[1] && > x<=fx[2]) retval <- 1 > else if (x>fx[2]) retval <- 0 > retval > } > > i(.4) > [1] 1 > > i(.3) > Error in i(0.3) : object 'retval' not > found > > When you convert your code for production use, look > at the ifelse and approx functions for doing this > sort of thing. They are much faster than > Vector(yourIntegrand). > > Also, when writing for help from R-help, please use > R syntax when defining variables, like > fx<-c(0.3,0.5,0.5,0.6) > and not some other language, like > fx<-[0.3,0.5,0.5,0.6] > If we can copy and paste your example into R instead > of translating the syntax by hand we will > be much more inclined to try to fix things up. > > Bill Dunlap > Spotfire, TIBCO Software > wdunlap tibco.com > > > -----Original Message----- > > From: r-help-boun...@r-project.org > > > [mailto:r-help-boun...@r-project.org] > On Behalf Of Julio Rojas > > Sent: Friday, December 18, 2009 7:02 AM > > To: r-help@r-project.org > > Subject: [R] Numerical Integration > > > > Dear @ll. I have to calculate numerical integrals for > > > triangular and trapezoidal figures. I know you can > calculate > > the exactly, but I want to do it this way to learn how > to > > proceed with more complicated shapes. The code I'm > using is > > the following: > > > > integrand<-function(x) { > > print(x) > > if(x<fx[1]) return(0) > > if(x>=fx[1] && x<fx[2]) > return((x-fx[1])/(fx[2]-fx[1])) > > if(x>fx[2] && x<=fx[3]) > return(1) > > if(x>fx[3] && x<=fx[4]) > return((x-fx[4])/(fx[3]-fx[4])) > > if(x>fx[4]) return(0) > > > > } > > > > fx<-data[i,j,] > > reltol<-1e-07 > > > integrate(Vectorize(integrand),0,1,rel.tol=reltol,subdivisions > > =200)$value > > > > It works for most cases, but then, I tried for the > triangle > > fx<-[0.3,0.5,0.5,0.6] and the following error was > presented: > > > > > > integrate(Vectorize(integrand),0,1,rel.tol=reltol,subdivisions=200) > > [1] 0.5 > > [1] 0.01304674 > > [1] 0.9869533 > > [1] 0.06746832 > > [1] 0.9325317 > > [1] 0.1602952 > > [1] 0.8397048 > > [1] 0.2833023 > > [1] 0.7166977 > > [1] 0.4255628 > > [1] 0.5744372 > > [1] 0.002171418 > > [1] 0.9978286 > > [1] 0.03492125 > > [1] 0.9650787 > > [1] 0.1095911 > > [1] 0.8904089 > > [1] 0.2186214 > > [1] 0.7813786 > > [1] 0.3528036 > > [1] 0.6471964 > > Error in integrate(Vectorize(integrand), 0, 1, rel.tol > = reltol, > > subdivisions = 200) : > > evaluation of function gave a result of wrong > type > > > > Does anybody know what happened? Thanks in advance!!! > > > > > > > > > ______________________________________________________________ > > ______________________ > > ¡Obtén la mejor experiencia en la web! > > Descarga gratis el nuevo Internet Explorer 8. > > http://downloads.yahoo.com/ieak8/?l=e1 > > > > ______________________________________________ > > 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. > > > ____________________________________________________________________________________ ¡Obtén la mejor experiencia en la web! Descarga gratis el nuevo Internet Explorer 8. http://downloads.yahoo.com/ieak8/?l=e1 ______________________________________________ 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.