R语言求根

时间:2024-07-04 18:03:50

求根是数值计算的一个基本问题,一般采用的都是迭代算法求解,主要有不动点迭代法、牛顿-拉富生算法、割线法和二分法。

  • 不动点迭代法

    所谓的不动点是指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
}
}
}