R语言求根

合集 - R语言(6)

1.R语言的精度和时间效率比较（简单版）2018-06-04

2.R语言求根2018-06-05

3.R语言最优化(多维）2018-06-144.R语言最优化(一维)2018-06-135.R语言RODBC数据库操作2018-06-086.R语言数值积分2018-06-07

收起

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

• 不动点迭代法

　　　　所谓的不动点是指x=f(x)的那些点，而所谓的不懂点迭代法是指将原方程化为x=f(x)形式之后，下一步所用的x值为这一步的f(x)，这样的话就可以一直逼近我们需　　　　　　　　　　　　         要的x，即方程的根，但是这种方法可能不会收敛到方程的根，随着初始值选定的大小，可能会有发散的情况，因此需要谨慎使用。

　　

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

###不动点迭代法 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))，因此这种方法需要有求导表达式，限制了使用范围。

　　　　

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

###牛顿迭代法 ###示例方程 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，我是用相似比解出来的。

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

###割线法 ###割线法公式：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之间，之后只要不断逼近根就好了。

?

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

###二分法 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     }   } }