R語言求根

求根是數值計算的一個基本問題，一般采用的都是迭代算法求解，主要有不動點迭代法、牛頓-拉富生算法、割線法和二分法。

• 不動點迭代法

　　　　所謂的不動點是指x=f(x)的那些點，而所謂的不懂點迭代法是指將原方程化為x=f(x)形式之后，下一步所用的x值為這一步的f(x)，這樣的話就可以一直逼近我們需　　　　　　　　　　　　         要的x，即方程的根，但是這種方法可能不會收斂到方程的根，隨着初始值選定的大小，可能會有發散的情況，因此需要謹慎使用。

　　

###不動點迭代法
func1 <- function(x){return(exp(exp(-x)))}
fixpoint <- function(func, x0, tol=1e-8, max.iter=1e4){
  ###求根的函數func
  ###初始值x0
  ###允許誤差范圍tol
  ###最大循環次數max.iter
  x.old <- x0
  x.new <- x0
  for(i in 1:max.iter){
    x.new <- func1(x.old)
    if(abs(x.new - x.old) < tol && i<max.iter){
      cat('the iter time is',i,'\n')
      return(format(x.new,digits = 9))
    }
    x.old <- x.new
  }
  cat('bad start num')
}

• 牛頓-拉富生算法

　　　　所謂的牛頓-拉富生算法其實就是課本里面說的牛頓迭代法，也不是一個難的程序，主要思想就是x(n+1)=x(n)-f(x(n))/f`(x(n))，因此這種方法需要有求導表達式，限制了使用范圍。

　　　　

###牛頓迭代法
###示例方程
func1 <- function(x){return(log(x)-exp(-x))}
func2 <- function(x){return(1/x+exp(-x))}
Newton <- function(func1, func2, x0, tol=1e-8, max.iter=1e4){
  ###牛頓迭代法//
  ###輸入方程式func1
  ###輸入func1的導函數func2
  ###初始值x0
  ###誤差范圍tol
  ###最大迭代次數
  x <- x0
  for(i in 1:max.iter){
    x <- x - func1(x)/func2(x)
    if(abs(x0-x) < tol){
      cat('iter num is:',i,'\n')
      return(x)
    }
    x0 <- x
  }
    cat('this is a not good start num or the function is bad!')
}

• 割線法

　　　　使用牛頓法的時候 ，需要先得到函數的導函數，但是很多時候都沒有導函數的解析表達式，因此可以考慮一種替代導函數的近似，即割線法，使用了一段直線作為函數某一段的近似，具體的數學表達式推導很簡單，在此略過，其主要假設是所求點的y值恰好為0，我是用相似比解出來的。

###割線法
###割線法公式：x(n+1)=[f(xn)x(n-1) - x(n)f(x(n-1))] / [f(xn) - f(x(n-1))]
###示例函數
func1 <- function(x){return(log(x)-exp(-x))}
gexian <- function(func,x0,x1,tol=1e-8,max.iter=1e4){
  ###割線法需要兩個初始值
  for(i in 1:max.iter){
    x <- (func(x1)*x0 - x1*func(x0)) / (func(x1) - func(x0))
    if(abs(x-x1) < tol){
      cat('iter num is:',i,'\n')
      return(x)
    }
    x0 <- x1
    x1 <- x
  }
  cat('this is a not good start num or the function is bad!')
}

•  二分法

　　　　上述的方法均有可能不收斂，並且對函數會有一些額外的要求，一種比較暴力的方法但是適用范圍更廣的解決辦法就是二分法，其思想是如果f(x)連續，那么就一定會存在f(x1)<0,f(x2)>0的情況，根必然在x1和x2之間，之后只要不斷逼近根就好了。

###二分法
func1 <- function(x){return(log(x)-exp(-x))}
half_M <- function(func,x0,x1,tol=1e-8,max.iter=1e4){
  ###初始x0,x1必須滿足f(x0)f(x1) < 0,x0 < x1
  for(i in 1:max.iter){
    cat(x0,"  ",x1,"\n")
    if(abs(x0-x1) < tol){
      cat('iter num:',i,"\n")
      return(x0)
    }
    temp_x <- (x0 + x1)/2
    if((func(temp_x)*func(x0)) < 0){
      x1 <- temp_x
    }
    else{
      x0 <- temp_x
    }
  }
}