ニュートン法で平方根の値を求める

きっかけ

電卓のルートってどうやって求めるんだっけ？　→　たぶんニュートン法だろうなぁ。ということで、Rで簡単に実装して計算してみました。

実装

ニュートン法

fnewton <- function(xn=1,p=2,a=2){
  return ( xn - ( xn^p - a )/( p * ((xn)^(p-1)) ) )
}

反復計算からグラフ作成

solve_with_fnewton <- function(a, p = 2, n, xn = 1 ){
  xs <- c(xn,xn+1:n)
  for(i in xs){
    xs[i+1] <- fnewton(xs[i],p,a)
  }
  last_xs <- xs[n]
  
  plot(xs,type="o")
  title(main=bquote(sqrt(.(a), .(p)) %~~% .(last_xs) ),sub = "")
  abline(h=a^(1/p) ,col="red")
  print( xs )
}

関数呼び出し

solve_with_fnewton(a = 2,n = 10,p = 2)

結果

感想

6回目の反復あたりでほぼ答えが電卓と一致して驚いた。