はじめに

今日で令和元年ゴールデンウィークも無事終了です。

前回に引き続きニュートン法で平方根を求めてみました。今回は、N = [1,1000] の範囲でニュートン法を用いて10回反復近似計算する時、真の値と近似値との残差はどんな感じになるのか、を調べてみます。

全てのNに対して初期値１から始めると、Nが大きくなればなるほど収束が遅れるような気配もありますが、まぁそこは部屋の片隅に置いておきましょう。

ニュートン法による平方根の近似

aの平方根を求めるには、初期値 $x_0$ を設定し、

$\displaystyle x_{n+1} \leftarrow x_n - \frac{ x_n^{2} - a }{ 2 x_n}$

を順次反復計算します。今回は $x_0 = 1$ とおきました。

例

a = 2, $x_0 = 1$ の時、

$\displaystyle x_{1} \leftarrow (1) - \frac{ ( 1 ) ^{2} - 2 }{ 2 \times ( 1 )} = 1.5$

$\displaystyle x_{2} \leftarrow (1.5) - \frac{ ( 1.5 ) ^{2} - 2 }{ 2 \times ( 1.5 )} = 1.416666667$

$\displaystyle x_{3} \leftarrow (1.416666667) - \frac{ ( 1.416666667 ) ^{2} - 2 }{ 2 \times ( 1.416666667 )} = 1.4142158687$

$\displaystyle x_{4} \leftarrow (1.4142158687) - \frac{ ( 1.4142158687 ) ^{2} - 2 }{ 2 \times ( 1.4142158687 )} = 1.414213563$

結論

初期値が1であっても、７回ほど反復計算しておけば問題なさそう。

コード

#rm(list=ls())

# f(x) = x^2 - a をニュートン法で求める更新式
fnewton <- function(xn=1,a=2){
  return ( xn - ( xn^2 - a )/( 2 * xn ) )
}

# 初期化
how_many_times_rep <- 10
N <- 1000
xs <- 1:N
root_xs <- sqrt(xs)
zeros <- rep(0,length(xs))

# 結果格納用行列
mat.x <- cbind( 
  matrix(c( xs,root_xs ),nrow = length(xs))
  ,matrix( rep(zeros, how_many_times_rep), nrow = length(zeros))
)

mat.x[,3] <- fnewton(xn = 1, a = mat.x[,1])

# ニュートン法にて更新
for(i in 4:(how_many_times_rep+2) ){
  mat.x[,i] = fnewton(xn = mat.x[,i-1], a = mat.x[,1])
}

#mat.dif <- matrix(0, nrow = nrow(mat.x), ncol = ncol(mat.x))
mat.dif <- mat.x

# 真の値 - 近似値
for(i in 3:(how_many_times_rep+2)){
  mat.dif[,i] <- mat.x[,2] - mat.x[,i]
}

# 結果確認
summary(mat.dif[,-c(1,2)])

par(family="Osaka")
boxplot(mat.dif[,-c(1,2)],
        col="lightpink",
        names=paste0(1:how_many_times_rep,"回目"),
        las=2,
        main=paste0("真の値と近似値との残差の箱ひげ図(N = 1 ~ ",N,")")
)

結果確認

要約統計量

summary の値を書き出すと、

V1               V2               V3               V4                V5         
 Min.   :-468.9   Min.   :-219.6   Min.   :-95.99   Min.   :-36.103   Min.   :-9.6226  
 1st Qu.:-348.2   1st Qu.:-161.4   1st Qu.:-69.00   1st Qu.:-24.697   1st Qu.:-5.8550  
 Median :-228.4   Median :-104.0   Median :-42.79   Median :-14.052   Median :-2.7105  
 Mean   :-229.7   Mean   :-105.3   Mean   :-44.04   Mean   :-15.158   Mean   :-3.3997  
 3rd Qu.:-110.0   3rd Qu.: -48.1   3rd Qu.:-18.09   3rd Qu.: -4.824   3rd Qu.:-0.5632  
 Max.   :   0.0   Max.   :   0.0   Max.   :  0.00   Max.   :  0.000   Max.   : 0.0000  

       V6                  V7                   V8                   V9            
 Min.   :-1.122493   Min.   :-1.924e-02   Min.   :-5.849e-06   Min.   :-5.400e-13  
 1st Qu.:-0.515571   1st Qu.:-4.763e-03   1st Qu.:-4.140e-07   1st Qu.:-3.553e-15  
 Median :-0.146451   Median :-4.762e-04   Median :-5.069e-09   Median : 0.000e+00  
 Mean   :-0.293577   Mean   :-3.268e-03   Mean   :-5.871e-07   Mean   :-3.005e-14  
 3rd Qu.:-0.009673   3rd Qu.:-2.953e-06   3rd Qu.: 0.000e+00   3rd Qu.: 0.000e+00  
 Max.   : 0.000000   Max.   : 0.000e+00   Max.   : 0.000e+00   Max.   : 3.553e-15  
      V10            
 Min.   :-3.553e-15  
 1st Qu.: 0.000e+00  
 Median : 0.000e+00  
 Mean   : 3.242e-17  
 3rd Qu.: 0.000e+00  
 Max.   : 3.553e-15

でした。

V6列(6回目の反復)で残差の絶対値の最小値がほぼ1
7回目で　　〃　　ほぼほぼ0

ですね。7回電卓ぽちぽちすれば平方根が得られそうです。

箱ひげ図

f:id:axjack:20190506230038p:plain

残差の範囲は反復回数が増す毎にどんどん小さくなっていきますね。５回目あたりで箱が見えなくなり、６回目あたりでヒゲもぺっちゃんこに見えます。

おわりに

計算方法覚えておけばルート機能要らない気もしないでも無いですね。脳内記憶容量を節約するならやはりルート付き電卓を買いましょう。

axjack's blog

### axjack is said to be an abbreviation for An eXistent JApanese Cool Klutz ###

初期値1から始めたニュートン法による平方根の近似値計算は、７回ほど反復計算しておけば大丈夫そうです。