Assignment 3

Course : STAT_6601

Assignment : 3

Submitted to : Prof John Angell

Submitted By : Devansh Pandit

NetID : vs5758

Instructions

Write two functions, root.bisect, and root.regula.falsi. The arguments to both functions are:

f, the function whose root will be found a and b, which are numbers such that a and b bracket a root of f. The function should check if a and b do bracket a root. If they do not, the function should stop with an appropriate message.

tol, a number which is used to determine how much accuracy is desired by the caller. This argument should have a default value.

You may have additional arguments, but you must explain what they are and their purpose.

The function should return the root found, the values of f evaluated at the root, and the number of iterations required to find the root.

The functions should have header comments to describe the arguments, and the output of the function, as well the purpose of the function. The rest of the code should be well-commented.

The function root.bisect will find the root using the bisection algorithm. The function root.regula.falsi will find the root using the regula.falsi algorithm. Do not use the R function uniroot or any other R function to find the roots. Your functions should not plot anything

##################################################
#
#    author: Devansh Pandit
#
#    function:    root.bisect
#
#    input arguments
#
#       f,           The function whose root we will calculate
#
#       x0,          A scalar
#       x1,          Another scalar. x0 and x1 are initial points
#                    for starting the algorithm
#      tol,          The algorithm ceases when the left and
#                    right points of the last brackting interval are
#                    within tol of each other.
#      max.iter,     The maximum number of iterations
#      .
#
#   returns:
#     a list containing $roots and $fvals as named list items.
#     Each of these is a vector. Naturally $roots contains the
#     root estimates, and $fvals contains the associated function
#     value at the estimated root.
#
#   
##################################################

##----------------------- Defining function f and function g
f <- function(x) { return(2*x*cos(2*x) - (x+1)^2) }
g<- function(x) { return( sin(x^(1/3))+ exp(-3*x) ) }


####----------------------start of bisection method

root.bisect <- function(f,x0,x1,tol = 1e-16,max.iter = 200)
{ # begin root.bisect
  
  if (!is.function(f))                   { stop('f is not a function.')}#f is function or not
  if (!is.vector(x0, mode = 'numeric'))  { stop('x0 is not numeric.')} #vector check
  if (!is.vector(x1, mode = 'numeric'))  { stop('x1 is not numeric.')} #vector check
  if ( length(x0) != 1 )                 { stop('x0 must have length equal to 1.')}
  if ( length(x1) != 1 )                 { stop('x1 must have length equal to 1.')}  
  if ( tol <  1.0e-17)                   { tol = 1.0e-17} #if tol less -17 make it -17
  if ( !is.numeric(max.iter) )           { stop( 'max.iter must be numeric.')}
  if ( max.iter < 1 )                    { stop( 'max.iter must be positive.')}
  
  
  iter.count = 0
  
  a=x0
  b=x1
  i=0
  tol=tol
  if (f(a)*f(b) > 0 ) { stop('f  function do not bracket root beween a and b.')}
  options(digits=20)
  while((abs(b-a) > tol) && (i < max.iter) )  ###((a+b)/2)!=a)
  { fa=f(a)
  fb=f(b)
  i=i+1
  m=0.5*(a+b)
  if ((m==a) || (m==b)){ print("root and any point matches") ; break} #if m is equal to a or b, means we get the root
  fm=f(m)    # calculating f of m
  if(f(m)*f(b) <0)  # checking if root can be found between values of these functions
  {a=m} # assigning new value to a if root exist between bracket intrval of m and b
  else
  {b=m  # assigning new value to b if root do not exist between bracket intrval of m and b
  }}
  
  return( list('roots' = m, 'f.vals' = f(m), 'iter.count' = i, 'tol'=tol))#returning values
}

# end of bisection method
  

  
  
####--------start of regulla falsi

root.regula.falsi <- function(f,x0,x1,tol = 1e-16,max.iter = 200)
{ # begin root.bisect
  
  if (!is.function(f))                   { stop('f is not a function.')}
  if (!is.vector(x0, mode = 'numeric'))  { stop('x0 is not numeric.')}
  if (!is.vector(x1, mode = 'numeric'))  { stop('x1 is not numeric.')}
  if ( length(x0) != 1 )                 { stop('x0 must have length equal to 1.')}
  if ( length(x1) != 1 )                 { stop('x1 must have length equal to 1.')}  
  if ( tol <  1.0e-17)                   { tol = 1.0e-17}
  if ( !is.numeric(max.iter) )           { stop( 'max.iter must be numeric.')}
  if ( max.iter < 1 )                    { stop( 'max.iter must be positive.')}
  
  
  iter.count = 0
  
  a=x0
  b=x1
  i=0
  tol=tol
  if (f(a)*f(b) > 0 ) { stop('f  function do not bracket root beween a and b.')}
  options(digits=20)
  while((abs(b-a) > tol) && (i < max.iter) )
  { fa=f(a)
  fb=f(b)
  i=i+1
  m=(((a*fb)-(b*fa))/(fb-fa))      ##### calculating X2 as m using formula for Line
  if (m==a|| (m==b)){ print("root and any point matches") ; break}
  fm=f(m) 
  if(f(m)*f(b) <0)
  {a=m}
  else
  {b=m
  }
  }
  
  return( list('roots' = m, 'f.vals' = f(m), 'iter.count' = i, 'tol'=tol ) )
}

# end of regula falsi




####--------start of regulla falsi modified

root.regula.falsi.modified <- function(f,x0,x1,tol = 1e-16,max.iter = 200)
{ # begin root.bisect
  
  if (!is.function(f))                   { stop('f is not a function.')}
  if (!is.vector(x0, mode = 'numeric'))  { stop('x0 is not numeric.')}
  if (!is.vector(x1, mode = 'numeric'))  { stop('x1 is not numeric.')}
  if ( length(x0) != 1 )                 { stop('x0 must have length equal to 1.')}
  if ( length(x1) != 1 )                 { stop('x1 must have length equal to 1.')}  
  if ( tol <  1.0e-17)                   { tol = 1.0e-17}
  if ( !is.numeric(max.iter) )           { stop( 'max.iter must be numeric.')}
  if ( max.iter < 1 )                    { stop( 'max.iter must be positive.')}
  
  
  iter.count = 0
  
  a=x0
  b=x1
  change.a=FALSE
  change.b=FALSE
  i=0
  fa=f(a)
  fb=f(b)/2
  tol=tol
  if (f(a)*f(b) > 0 ) { stop('f  function do not bracket root beween a and b.')}
  options(digits=15)
  while((abs(b-a) > tol) && (i < max.iter) )  ###((a+b)/2)!=a)
  { 
    #print(i)  
    i=i+1
    m=(((a*fb)-(b*fa))/(fb-fa)) #(((a*fb)-(b*fa))/(fb-fa))      ##### calculating X2 as m
    if (m==a || (m==b)){break}
    
    fm=f(m)
    
    #testing
    if(f(m)*f(b) <0)
    {a=m
    if (change.a){fb = fb/2 ; change.a=TRUE ;change.b=FALSE}  # making function as half
    #making funtion half is the base for modified regulla falsi method
    fa=f(a)
    }
    else
    {b=m
    #fa = fa/2
    if (change.b){fa = fa/2 ; change.b=TRUE ;change.a=FALSE}
    fb=f(b)
    }
  }
  
  
  return( list('roots' = m, 'f.vals' = f(m), 'iter.count' = i, 'tol'=tol ) )
}

# end of regula falsi Modified

Question 1 (1) Find the root of the function

f <- function(x) { return(2xcos(2*x) - (x+1)^2) } for the following intervals: first, a = -3, b = -2. Second, a = -2, b = 4 Report the value of tol used, the roots found, the values of f evaluated at the roots found, and the number of iterations. Run both functions.

print("root information from Bisection method")

## [1] "root information from Bisection method"

root.bisect(f,-3,-2,tol=1e-15,max.iter=70)

## $roots
## [1] -2.1913080117972461
## 
## $f.vals
## [1] 5.773159728050814e-15
## 
## $iter.count
## [1] 50
## 
## $tol
## [1] 1.0000000000000001e-15

root.bisect(f,-2,4,tol=1e-15,max.iter=70)

## $roots
## [1] -0.7981599614057957
## 
## $f.vals
## [1] -7.5633943552588789e-16
## 
## $iter.count
## [1] 53
## 
## $tol
## [1] 1.0000000000000001e-15

try(root.bisect(f,0,5,tol=1e-15,max.iter=70))

print("root information from Regulla Falsi method")

## [1] "root information from Regulla Falsi method"

root.regula.falsi(f,-3,-2,tol=1e-15,max.iter=70)

## $roots
## [1] -2.191308011797247
## 
## $f.vals
## [1] -3.1086244689504383e-15
## 
## $iter.count
## [1] 20
## 
## $tol
## [1] 1.0000000000000001e-15

root.regula.falsi(f,-2,4,tol=1e-15,max.iter=70)

## [1] "root and any point matches"

## $roots
## [1] -0.79815996140579593
## 
## $f.vals
## [1] 4.8572257327350599e-17
## 
## $iter.count
## [1] 35
## 
## $tol
## [1] 1.0000000000000001e-15

print("root information from Modified Regulla Falsi method")

## [1] "root information from Modified Regulla Falsi method"

root.regula.falsi.modified(f,-3,-2,tol=1e-15,max.iter=70)

## $roots
## [1] -2.19130801179725
## 
## $f.vals
## [1] 1.33226762955019e-15
## 
## $iter.count
## [1] 21
## 
## $tol
## [1] 1e-15

root.regula.falsi.modified(f,-2,-4,tol=1e-15,max.iter=70)

## $roots
## [1] -2.19130801179725
## 
## $f.vals
## [1] 1.33226762955019e-15
## 
## $iter.count
## [1] 23
## 
## $tol
## [1] 1e-15

Question 2 . Now call the two functions to find a root of f in the interval a = 0, b = 5. Report what happens.

Note: I have taken answer from console, as it is giving error in R Markdown. We are not getting root from bracket interval 0 and 5 and as per error check, we are getting error.

##root.bisect(f,-3,-2,tol=1e-15,max.iter=70)


##root.regula.falsi(f,-3,-2,tol=1e-15,max.iter=70)

Question 3 Find the first three positive roots (that is, the three smallest positive roots) of the function

g<- function(x) { return( sin(x^(1/3))+ exp(-3*x) ) } Use both functions to solve this problem. Report the bracketing intervals used as well as the information reported for problem (1). The bracketing intervals should have a length of at least 4. You may plot the function g to find appropriate bracketing intervals, or you may simply evaluate it at random points until you find a bracketing interval.

g<- function(x) { return( sin(x^(1/3))+ exp(-3*x) ) }
a=1
b=5
for ( i in 1:4)  # can put any number instead of 1:4 , like 1:10 will make 10 roots
{
while ( g(a)*g(b)> 0)  ## while loop will do looping until it get bracket interval for root
{
  a=a+1
  b=b+1
  
}
cat("\n")  
print(paste0("Bracketing interval : ",i))
print (a)
print(b)
cat("\n") 
cat("\n") 
print(paste0("smallest positive roots : ",i))
cat("\n") 
print(paste0("smallest positive roots from bisection algorithm : ",i))
print(root.bisect(g,a,b,tol=1e-15,max.iter=70))

cat("\n") 
cat("\n") 
print(paste0("smallest positive roots from regulla falsi algorithm : ",i))
print(root.regula.falsi(g,a,b,tol=1e-15,max.iter=70))

cat("\n") 
cat("\n") 
print(paste0("smallest positive roots from modified regulla falsi algorithm : ",i))
print(root.regula.falsi.modified(g,a,b,tol=1e-15,max.iter=70))

a=a+4
b=b+4
}

## 
## [1] "Bracketing interval : 1"
## [1] 28
## [1] 32
## 
## 
## [1] "smallest positive roots : 1"
## 
## [1] "smallest positive roots from bisection algorithm : 1"
## [1] "root and any point matches"
## $roots
## [1] 31.006276680299827
## 
## $f.vals
## [1] 1.2246063538223773e-16
## 
## $iter.count
## [1] 51
## 
## $tol
## [1] 1.0000000000000001e-15
## 
## 
## 
## [1] "smallest positive roots from regulla falsi algorithm : 1"
## [1] "root and any point matches"
## $roots
## [1] 31.006276680299823
## 
## $f.vals
## [1] 1.2246063538223773e-16
## 
## $iter.count
## [1] 11
## 
## $tol
## [1] 1.0000000000000001e-15
## 
## 
## 
## [1] "smallest positive roots from modified regulla falsi algorithm : 1"
## $roots
## [1] 31.0062766802998
## 
## $f.vals
## [1] 1.22460635382238e-16
## 
## $iter.count
## [1] 13
## 
## $tol
## [1] 1e-15
## 
## 
## [1] "Bracketing interval : 2"
## [1] 245
## [1] 249
## 
## 
## [1] "smallest positive roots : 2"
## 
## [1] "smallest positive roots from bisection algorithm : 2"
## [1] "root and any point matches"
## $roots
## [1] 248.05021344239867
## 
## $f.vals
## [1] 6.4325714893564978e-16
## 
## $iter.count
## [1] 48
## 
## $tol
## [1] 1.0000000000000001e-15
## 
## 
## 
## [1] "smallest positive roots from regulla falsi algorithm : 2"
## [1] "root and any point matches"
## $roots
## [1] 248.05021344239861
## 
## $f.vals
## [1] -2.4492127076447545e-16
## 
## $iter.count
## [1] 7
## 
## $tol
## [1] 1.0000000000000001e-15
## 
## 
## 
## [1] "smallest positive roots from modified regulla falsi algorithm : 2"
## $roots
## [1] 248.050213442399
## 
## $f.vals
## [1] 6.4325714893565e-16
## 
## $iter.count
## [1] 9
## 
## $tol
## [1] 1e-15
## 
## 
## [1] "Bracketing interval : 3"
## [1] 834
## [1] 838
## 
## 
## [1] "smallest positive roots : 3"
## 
## [1] "smallest positive roots from bisection algorithm : 3"
## [1] "root and any point matches"
## $roots
## [1] 837.16947036809552
## 
## $f.vals
## [1] 3.6738190614671318e-16
## 
## $iter.count
## [1] 46
## 
## $tol
## [1] 1.0000000000000001e-15
## 
## 
## 
## [1] "smallest positive roots from regulla falsi algorithm : 3"
## [1] "root and any point matches"
## $roots
## [1] 837.16947036809563
## 
## $f.vals
## [1] -1.4089749332535373e-15
## 
## $iter.count
## [1] 7
## 
## $tol
## [1] 1.0000000000000001e-15
## 
## 
## 
## [1] "smallest positive roots from modified regulla falsi algorithm : 3"
## $roots
## [1] 837.169470368096
## 
## $f.vals
## [1] -1.40897493325354e-15
## 
## $iter.count
## [1] 8
## 
## $tol
## [1] 1e-15
## 
## 
## [1] "Bracketing interval : 4"
## [1] 1981
## [1] 1985
## 
## 
## [1] "smallest positive roots : 4"
## 
## [1] "smallest positive roots from bisection algorithm : 4"
## [1] "root and any point matches"
## $roots
## [1] 1984.4017075391894
## 
## $f.vals
## [1] -4.898425415289509e-16
## 
## $iter.count
## [1] 45
## 
## $tol
## [1] 1.0000000000000001e-15
## 
## 
## 
## [1] "smallest positive roots from regulla falsi algorithm : 4"
## [1] "root and any point matches"
## $roots
## [1] 1984.4017075391894
## 
## $f.vals
## [1] -4.898425415289509e-16
## 
## $iter.count
## [1] 6
## 
## $tol
## [1] 1.0000000000000001e-15
## 
## 
## 
## [1] "smallest positive roots from modified regulla falsi algorithm : 4"
## $roots
## [1] 1984.40170753919
## 
## $f.vals
## [1] -4.89842541528951e-16
## 
## $iter.count
## [1] 7
## 
## $tol
## [1] 1e-15

########################## funtion to find bracketing interval