Bisect(二分法)

cubic <- function (x, a=0, b=0, c=0, d=0) a*x^3 + b*x^2 + c*x + d
a <- 1 ; b <- -2 ; c <- 3 ; d <- -4
result <- uniroot(cubic, c(-5,5), a = a, b = b,c = c,d = d, tol = 0.000001)
result

## $root
## [1] 1.650629
## 
## $f.root
## [1] 5.33745e-09
## 
## $iter
## [1] 8
## 
## $init.it
## [1] NA
## 
## $estim.prec
## [1] 5e-07

# $`root`
# [1] 1.650629
# 
# $f.root
# [1] 5.33745e-09
# 
# $iter
# [1] 8
# 
# $init.it
# [1] NA
# 
# $estim.prec
# [1] 5e-07

bisect <- function(f, a, b, eps = 1e-5, ...){
  if (!is.function(f)){
    list(fail="f is not a function!")
  }
  else if(f(a, ...) * f(b, ...) > 0)
    list(fail="finding root is fail!")
  else{
    repeat{
    if (abs(b - a) < eps) break
    x <- (a + b)/2
    if (f(a, ...) * f(x, ...) < 0) b <- x else a <- x
    }
    list(root = (a+b)/2, fun = f(x, ...))
}
}
bisect(cubic, -5, 5, eps = 1e-6, 1, -2 , 3, -4)

## $root
## [1] 1.650629
## 
## $fun
## [1] 2.048755e-06

# $`root`
# [1] 1.650629
# 
# $fun
# [1] 2.048755e-06

Newton(牛顿法)

Intro

If we want to solve the system of nonlinear equations，then \[ x^{k+1} = x^{(k)}-[J(x^{(k)})]^{-1}f(x^),k = 0,1,\cdots \]

And Jacobi matrix:

\[\begin{equation} \textbf{J}(x) = \begin{bmatrix} \dfrac{\partial f_1}{\partial x_1} & \dfrac{\partial f_1}{\partial x_2} & \cdots &\dfrac{\partial f_1}{\partial x_n} \\ \dfrac{\partial f_2}{\partial x_1} & \dfrac{\partial f_2}{\partial x_2} &\cdots& \dfrac{\partial f_2}{\partial x_n} \\ \vdots& \vdots& \ddots & \vdots \\ \dfrac{\partial f_m}{\partial x_1} &\dfrac{\partial f_m}{\partial x_2} & \cdots &\dfrac{\partial f_m}{\partial x_n} \\ \end{bmatrix} \end{equation}\]

If we want to solve this equation (start = c(0,1), ep <= 1e-5):

\[\begin{equation} \left\{ \begin{aligned} \ x_1^2 + x_2^2 - 5 = 0\\ \ (x_1 + 1)x_2 - (3x_2 + 1) = 0\\ \end{aligned} \right. \end{equation}\]

Newtons<-function (fun, x, ep = 1e-5, it_max = 100){
  index <- 0; k <- 1
  while (k <= it_max){
    x1 <- x; obj <- fun(x);
    x <- x - solve(obj$J, obj$f);
    norm <- sqrt((x - x1) %*% (x - x1))
    if (norm < ep){
      success <- 1; break
  }
    k <- k + 1
  }
  obj <- fun(x);
  list(root = x,  val = obj$f, iteration = k, success = success)
}

funs <- function(x){
  f <- c(x[1]^2 + x[2]^2 - 5, (x[1] + 1) * x[2] - (3 * x[1] + 1))
  J <- matrix(
    c(2 * x[1], 2 * x[2], x[2] - 3, x[1] + 1),
    nrow = 2, byrow = T)
  list(f = f, J = J)
}

Newtons(funs, c(0,1))

## $root
## [1] 1 2
## 
## $val
## [1] 1.598721e-14 6.217249e-15
## 
## $iteration
## [1] 6
## 
## $success
## [1] 1

Differences(差分)

diff(x, lag = 1, differences = 1, …) lag(x, k = 1, …)

diff(1:10)

## [1] 1 1 1 1 1 1 1 1 1

#  1 1 1 1 1 1 1 1 1
diff(1:10, lag=2,differences=1)

## [1] 2 2 2 2 2 2 2 2

# 2 2 2 2 2 2 2 2
diff(1:10, 2, 2)

## [1] 0 0 0 0 0 0

# 0 0 0 0 0 0
lag(1:10,5)

##  [1]  1  2  3  4  5  6  7  8  9 10
## attr(,"tsp")
## [1] -4  5  1

#  [1]  1  2  3  4  5  6  7  8  9 10
# attr(,"tsp")
# [1] -4  5  1
lagdf <- function(x, k) {
    c(rep(NA, k), x)[1 : length(x)] 
}

Integral(积分)

Intro

Newton-Leibniz Formula:

\[ \int_{a}^{b} f(x)dx = F(b) - F(a) \]

梯形公式(n=1):

\[ T= (b - a)[\frac {1}{2} f(a) + \frac {1}{2} f(b)] \]

辛普森公式(n=2, Simpson’s rule):

\[ S=(b-a)[\frac {1}{6} f(a) + \frac {4}{6} f(\frac{a+b}{2}) + \frac {1}{6} f(b) ] \] \(Divide [a,b] to\, n\, fractions, step\, h = \frac{b-a}{n}, then\, x_i = a + ih,i = 1,2,\cdots,x_n:\)

复化梯形公式: \[ \int_{a}^{b} f(x)dx \approx \sum_{i=1}^{n}\int_{x_i}^{x_{i+1}}dx=T_n= \frac{h}{2}[f(a) + 2\sum_{i=1}^{n-1}f(a+ih) + f(b)] \] 复化辛普森公式: \[ S_n = \frac{h}{6}\sum_{i=0}^{n-1}[f(x_i) + 4f(x_{i + \frac{1}{2}}) + f(x_{i+1}))] \]

复化梯形公式

area <- function(f, a, b, eps = 1.0e-06, lim = 10) {
  fun1 <- function(f, a, b, fa, fb, a0, eps, lim, fun) {
  d <- (a + b) / 2
  h <- (b - a) / 4
  fd <- f(d)
  a1 <- h * (fa + fd)
  a2 <- h * (fb + fd)
  if (abs(a0 - a1 - a2) < eps || lim == 0)
    return(a1 + a2)
  else {
    return(
      fun(f, a, d, fa, fd, a1, eps, lim - 1, fun) + 
      fun(f, d, b, fd, fb, a2, eps, lim - 1, fun)
      )
  }
  }
  fa <- f(a)
  fb <- f(b)
  a0 <- ((fa + fb) * (b - a)) / 2
  fun1(f, a, b, fa, fb, a0, eps, lim, fun1)
  }

f <- function(x) 1 / x^2
  
quad <- area(f, 1, 10)
quad