Question no. 1

func <- function(x) {
  x - 25^1/3
}

curve(func, xlim=c(6,10), col='blue', lwd=1.5, lty=2)
abline(h=0)
abline(v=0)

Lets define bisection method as:

bisection <- function(f, a, b, n = 1000, tol = 1e-7) {
  if (f(a)*f(b) > 0) {
    stop('signs of f(a) and f(b) does not differ')
  }
  
  for (i in 1:n) {
    c <- (a + b) / 2
    
    if ((f(c) == 0) || ((b - a) / 2) < tol) 
      return(c)
    
    
    if(f(c)*f(a) < 0)
      b <- c
    
    else
      a <- c
  }
  print('Too many iterations')
}

Lets choose a = 8, b = 9, since the root is in this boundary. And also 1e-4 as tolerance as well.

bisection(func, a = 8, b = 9, tol = 1e-4)
## [1] 8.333313

Question no. 2

I defined the recursive function as like this:

a <- runif(1, min=.Machine$double.eps, max = 1e50)
a
## [1] 4.390204e+49
x0 <- runif(1, min=0, max=2/a)
x0
## [1] 4.52437e-51
func <- function(n){
  if (n == 0) 
    return(x0)
  
  else{
    xnm1 = func(n-1)
    return(xnm1*(2-a*xnm1))
  }
}

To show that the function has a limit, we should iterate throw different values of n starting from 1 to some larger numbers.

res <- c()

for (n in 1:20) {
  res[n] <- func(n)  
}

plot(res, type = "l")

From the plot above we can deduce that the function \(func\) rapidly converges to an specific number.
and we can get that number just by accessing the last element of vector \(res\): 2.2777983^{-50}

Question no. 9

A <- runif(1, min=.Machine$double.eps, max = 1e3)
x1 <- runif(1, min=.Machine$double.eps, max = 1e3)

func <- function(n){
  if (n == 1) 
    return(x1)
  
  else{
    xnm1 = func(n-1)
    return(0.5*xnm1 + A/(2*xnm1))
  }
}

To show that the function has a limit, we should iterate throw different values of n starting from 1 to some larger numbers.

res <- c()

for (n in 1:20) {
  res[n] <- func(n)  
}

plot(res, type = "l")
abline(h=sqrt(A), lwd=1.5, lty=2, col="red")

We can calculate the residual by the following code:

sqrt(A) - res[length(res)]
## [1] 0

From observance the function as convergence to \(\sqrt A\).

Question no. 8

First lets define equation of \(ln(x) = x^2 - 1\) as a function of \(f(x) = ln(x) - x^2 + 1\).

func <- function(x){
  return(log2(x) - x^2 +1)
}

curve(func, from = 0.1, to=6)
abline(h = 0, lwd=1.5, lty=2, col="red")

A little bit of zoom on roots…

curve(func, from = 0.4, to=1.4)
abline(h = 0, lwd=1.5, lty=2, col="red")
abline(v = 0.45, lwd=1.5, lty=2, col="blue")

It seems that 0.45 is not the root of the function!
But changing the log function’s base from \(e\) to \(10\) makes scene!

func <- function(x){
  return(log(x) - x^2 +1)
}

curve(func, from = 0.4, to=1.4)
abline(h = 0, lwd=1.5, lty=2, col="red")
abline(v = 0.45, lwd=1.5, lty=2, col="blue")

Lets define the iterative function of \(func\) to analysis fixed-point method.

ifunc <- function(x){
  if(x < 0){
    warning("Out of range of func.")
    return(NA)
  }  
  return(sqrt(1+log(x)))
}

Then we declare fixed-point method’s function as bellow:

fixedPoint <- function(x0, func, eps, color, holdPlot=F){
  
  if(holdPlot == F){
    curve(ifunc, 0.3, 1.2)
    lines(x=c(-10,10), y=c(-10,10), col="blue")
  }
  lines(x=c(x0, x0), y=c(0, ifunc(x0)), lty=2, col=color)
  
  prevx <- x0
  for (r in 1:50) {
    x <- ifunc(prevx)
    if(is.na(x)){
      break
    }
    
    lines(x=c(prevx, x), y=c(x, x), lty=2, col=color)
    lines(x=c(x, x), y=c(x, ifunc(x)), lty=2, col=color)
    
    
    if(abs(x - prevx) < eps){
      message(x)
      message("iterations spent: ", r)
      break
    }
    else
      prevx <- x
  }
}

Let’s test the functionality by passing the iterative function of \(ifunc\) to \(fixedPoint\) function, initializing variable \(x0\) as 1.1 and \(eps\) as 1e-5.

x0 <- 1.1
eps <- 1e-5

fixedPoint(x0, ifunc, eps, "red")
## Warning in if (x < 0) {: the condition has length > 1 and only the first element
## will be used
## Warning in sqrt(1 + log(x)): NaNs produced
## 1.00000531576379
## iterations spent: 14

We see the result converges to 1 after 14 iterations. Lets test by other initial values for \(x0\) and plot together.

colors <- c("red", "brown", "green3", "orange", "purple")
x0 <- c(0.44, 0.5, 0.6, 0.75, 1.2)

for (test in 1:length(x0)) {
  suppressWarnings(fixedPoint(x0[test], ifunc, eps, colors[test], ifelse(test == 1, F, T)))
}
## 0.999993478450699
## iterations spent: 19
## 0.999991261634415
## iterations spent: 17
## 0.99999373972343
## iterations spent: 16
## 1.00000943891634
## iterations spent: 14

As we saw in the above plot, no matters how much we get closer to the left-side-root, we always get converged to the right-side-root. Values smaller than 0.45 diverges.

Conclusion: left-side-root can not be calculated by fixed-point method.

Question no. 10

Let’s define the inside of integral as \(integrand(x)\)

integrand <- function(x){
  return(exp(-1*(x^2)))
}

curve(integrand, -5, 5)

Next we need a method to get the integrals answer. Paper in bellow represents an useful method:
http://homepages.math.uic.edu/~jyang06/stat522/handouts/handout8.pdf

by using the following command we can calculate the integral:

integrate(integrand, lower = 0, upper = Inf)
## 0.8862269 with absolute error < 2.2e-06

The absolute value can be accessed by ‘$value’.

Since we the calculation method is known, lets create a vector to store integral’s answers by some steps in an specific boundaries.

xs <- seq(from=-4, to=4, by=0.1)
ys <- numeric(length(xs))
i  <- 0
for (dx in xs) {
  ys[i] <- integrate(integrand, lower = 0, upper = dx)$value
  i <- i + 1
}
xs[length(xs)] <- NA
ys[length(ys)] <- NA
plot(xs, ys)
lines(xs, ys)

Then we need to multiply the function by 10 and subtract 1 to get the final function. The compressed construction would be:

func <- function(x){
  res <- integrate(function(x){
    exp(-1*(x^2))
  }, lower = 0, upper = x)$value
  return(10*res-1)
}

xs <- seq(from=-5, to=5, by=0.1)
ys <- numeric(length(xs))
i  <- 0
for (dx in xs) {
  ys[i] <- func(dx)
  i <- i + 1
}
xs[length(xs)] <- NA
ys[length(ys)] <- NA
plot(xs, ys)
lines(xs, ys)
abline(h = 0, lwd=1.5, lty=2, col="red")

I actually don’t now how to get the derivative of \(func\)!! So lets numerically deviate and plot it.

dys <- c(NA, ys[2:length(ys)] - ys[1:(length(ys)-1)])
plot(xs, dys)
lines(xs, dys)

I found this so close to the first \(integrand()\) function. Lets modify and plot them together.

dys <- c(NA, ys[2:length(ys)] - ys[1:(length(ys)-1)])/c(NA, xs[2:length(xs)] - xs[1:(length(xs)-1)])
plot(xs, dys)
lines(xs, dys)

derivate <- function(x){
  return(exp(-1*(x^2))*10)
}

fys <- c()
i <- 1
for (dx in xs) {
  fys[i] <- derivate(dx)
  i <- i + 1
}
lines(xs, fys, col="red")

Lets assume \(derivate()\) is the first deviate of the function \(func()\)