dimanche 15 novembre 2020

vectorized IF for bisection method in R

i m a beginner in R and i'm spending hours and hours trying to debug a function with IFELSE statement. I know it could be kind of obvious but i can't find a solution. Here' s my problem.

Using the same funcion with just an unknown factor everything works. here's my function:

tir <- function(i){
 return((-45000 * 1-( 1 + i)^-3)/i + (1600000 * (1+i)^-3))
 }
tir(0.2)

> tir(0.2)
[1] 700923

Then i solve it with the bisection method to find the root.

bisec5 <- function(fun, a, b, tol=0.00001) {
  if (abs(fun(a)) < tol )
    return(a)
  if (abs(fun(b)) < tol )
    return(b)
  x <- (a + b) / 2
  if( abs(fun(x)) < tol )
    return(x)
  if (fun(a)*fun(x) < 0)
    return( Recall( fun, a, x, tol) )
  else
    return( Recall( fun, x, b, tol) )
}

bisec5(fun = tir, a= 50, b= 0.00000001)

[1] 4.380243

My issues begin when i'm tryng to solve the same equation for a vector. here the first function when the exponent of (1+i) is a number included between 2 and 28

tir <- function(i){
    n<- 2:28
    return((-45000 * 1-( 1 + i)^-n)/i + (1600000 * (1+i)^-n))
}
tir(0.2)
[1]  886107.639  700923.032  546602.527  418002.106  310835.088  221529.240  147107.700   85089.750
 [9]   33408.125   -9659.896  -45549.913  -75458.261 -100381.884 -121151.570 -138459.642 -152883.035
[17] -164902.529 -174918.774 -183265.645 -190221.371 -196017.809 -200848.174 -204873.479 -208227.899
[25] -211023.249 -213352.707 -215293.923

Finally, how could i find the 28 roots of this equation applying the same method i used for 1? Using the exactly the same function i just find the first root. therefore, I've tried using IFELSE, since the solution im looking for will be a vector, but i keep find errors. here's one example of my many tryings.

bisec5 <- function(fun, a, b, tol=0.00001) {
  ifelse (abs(fun(a)) < tol, a, "0" )
    
  ifelse (abs(fun(b)) < tol, print(b), next)
   
  x <- (a + b) / 2
  ifelse ( abs(fun(x)) < tol, print(x), next )
   
  ifelse (fun(a)*fun(x) < 0,return( Recall( fun, a, x, tol)), return( Recall( fun, x, b, tol)))
}
bisec5(fun = tir, a= 50, b =0.000001)
Error in ifelse(abs(fun(b)) < tol, print(b), next) : 
  no loop for break/next, jumping to top leve

Hope i coluld explain my problem as clear as possible, be comprehensive to any lack or mistake of terminology i made, and thanks in advance.

Aucun commentaire:

Enregistrer un commentaire