在R中可以使用uniroot函數來求解一元方程。
R中uniroot函數的語法形式如下:
uniroot(f, interval, ...,
lower = min(interval), upper = max(interval),
f.lower = f(lower, ...), f.upper = f(upper, ...),
extendInt = c("no", "yes", "downX", "upX"), check.conv = FALSE,
tol = .Machine$double.eps^0.25, maxiter = 1000, trace = 0)
主要參數的含義如下:
(1)f 要求解的方程;
(2)interval 用向量表示的區間;
(3)lower,upper 用于搜索解的下限和上限;
(4)tol 求解的精度;
(5)maxiter 最大迭代次數
該函數的結果是一個列表,包括4個部分:求解的根root和在該點的函數值f.root;迭代次數iter和求解方程的近似估計的精度estim.prec。
下面舉幾個例子來說明一下具體使用方法:
(1)求解一元一次方程
求解形如ax+b=0的根。如求解 3x + 2 = 0的根:
f <- function(x, a, b) { return (a * x + b) }
root <- uniroot(f, c(-5, 0), a = 3, b = 2, tol=0.0001)
root$root #輸出根
root$f.root #輸出x值為求解的根時的函數值
以上運行結果如下圖所示:
上面的程序也可以寫成下面的形式:
root<-uniroot(function(x,a,b) a*x+b,a=3,b=2,lower=-5,upper=0,tol=0.0001)
結果是一樣的。
(2)求解一元二次方程
求解形如f(x)=ax2 + bx + c的方程的解。如求出方程x2 - x - 6 = 0的根。
f <- function(x, a, b, c) { return (a * x ^ 2 + b * x + c)}
re<-uniroot(f,c(-3,0),a=1,b=-1,c=-6,tol=1e-9)
re$root #輸出根
re$f.root #當方程取值根時的方程結果
re$iter #求解的迭代次數
運算結果如下圖所示:
從上圖可以看出該方程的一個根為-2,改變它的運算區間可以求出另外一個根:
re<-uniroot(f,c(0,3),a=1,b=-1,c=-6,tol=1e-9)
re$root
這里求得的另外一個根為3.
由于uniroot()函數每次只能計算一個根,而且要求輸入的區間端點值必須是正負號相反的。如果我們直接輸入(-3, 4)這個區間,那么uniroot()函數會出現錯誤:
Error in uniroot(f, c(-3, 4), a = 1, b = -1, c = -6, tol = 1e-09) :
位于極點邊的f()值之正負號不相反
即要求f(lower) * f(upper) <=0。
那么在實際中如何確定根的大致范圍呢,一個做法就是通過繪圖來觀察,如本例中可以使用西面的代碼來繪圖:
f<-function(x,a,b,c){return (a*x^2+b*x+c)}
X<-seq(-6,6,by=0.1)
Y<-f(X,a=1,b=-1,c=-6)
plot(X,Y,type="l",col=2) #繪制函數圖形
abline(h=0,v=0,lty=2,col="gray") #繪制出橫軸和縱軸
繪圖結果如下圖所示:
從圖中可以看出,方程的一個根位于-3和0之間,另外一個根位于2和4之間,這樣我們取區間(-3,0)和(0,4)即可。同時也可以看出這兩個區間的兩個端點的函數值得符號是相反的。
(3)求解一元三次方程
求解形如ax3+bx2+cx+d=0的根。如求x3-x-3=0的方程根。
程序如下:
f<-function(x,a,b,c,d) {return(a*x^3+b*x^2+c*x+d)}
a<-1
b<-0
c<- -1
d<- -3
X<-seq(-3,4,by=0.01)
Y<-f(X,a=a,b=b,c=c,d=d)
plot(X,Y,type="l",col=3)
abline(h=0,v=0,col="gray")
re<-uniroot(f,c(1,2.5),a=a,b=b,c=c,d=d,tol=1e-10)
re$root
re$f.root
re$iter
求解結果如下:
繪制的曲線如下圖所示:
從圖中可以看出方程的根在區間(1,3)之間,且求得的根為1.6717
新聞熱點
疑難解答