Rで普通の数値解析その１
塩には必要だが。
あっ, どーも僕です。
来年から
おかげさまで、来年からの食い扶持が決まりました。
特に@beroberoさんは見ず知らずの私に声をかけてくださりありがとうございました。励みになりました。
内定先ではデータ解析とかとはあまり関係ない仕事かもですが、分野としては志望通りで一安心です。
Rで普通の数値解析
近頃というか、前々からのですが、これからのためには普通の数値解析をちゃんと実装したことほうがいいのでわと感じていました。
そこで、今回はRを使ってですが普通の数値解析を勉強しました。今週を使ってコツコツやりました。ただ、ここで実装したのはすべてRで提供されていますのであくまで自分の訓練用です。
それなら、せめてRcpp*1使えって話ですがRcppの教科書は一回しか読んでなくてまだ使いこなせていのでやめました。
教科書はこちら。
- 作者: 金谷健一
- 出版社/メーカー: 共立出版
- 発売日: 2010/10/23
- メディア: 単行本
- クリック: 10回
- この商品を含むブログを見る
初読は一年くらい前ですが、そのときはただ読んだだけでした。内容のレベルはまえがきによると情報系の学部生向けです。
ちなみに数値解析をRでやるならpracma*2がオススメです。
ってことで、以下には各章ごとに実装したRのプログラムをいくつか載せます。気になる項目があったら見てみてください。一応、全てちゃんと動作しています。ただ、ピポットの選択は自信ないです。（教科書に例題がなかった）また、実装したの次のものです。ところどころ、章をまたいで依存してるかもです。
あと、誤差のことは説明するの大変なので、本買ってください。
第一章 数値と誤差
- 十進法を二進法に変換
- 二進法を十進法に変換
第五章 反復法と収束
- ヤコビ法
- ガウスサイデル（コツコツ）
- ガウスサイデル（行列使いつつ）
追記
スクリプトの最新版はこちらから。
Rで普通の数値解析
第一章 数値と誤差
十進法を二進法に変換
# *********************************************************************************************** # 10 → 2 my10to2 <- function (x) { x <- abs (x) z <- floor (x) # 整数部 d <- x - z # 小数部 ansZ <- numeric (1e3) ansD <- numeric (1e3) # 整数部の処理 res <- 1 n <- 1 # index if (z != 0) { while (res != 0) { res <- z %/% 2 m <- z %% 2 z <- res ansZ [n] <- m n <- n + 1 } } else { ansZ <- 0 } ansZ <- paste(ansZ[(n-1):1], collapse = "") # a_0から埋めていったので # 小数部の処理(小数点以下10桁まで) res <- 1 n <- 1 if (d != 0) { while (d != 0) { tmp <- 2 * d res <- floor (tmp) d <- tmp - res ansD[n] <- res n <- n + 1 if (n > 10) break } } else { ansD <- 0 } ansD <- paste (ansD[1:(n-1)], collapse = "") # 結果をまとめて出力 ans <- paste (ansZ, ansD, sep = ".") return (ans) } my10to2 (417.762) # 100001011.1100001100 my10to2 (4) # 100
二進法を十進法に変換
# 2 → 10 my2to10 <- function (x, digits = 4) { x <- abs (x) z <- floor (x) d <- x - z ansZ <- numeric(1) ansD <- numeric(1) # 整数部の処理 z <- as.character (z) z <- as.numeric(strsplit (z, "")[[1]]) ansZ[1] <- z[1] for (i in seq_len (length (z) - 1)) { ansZ <- crossprod (c(2, 1), z[c (i, i+1)]) z[i + 1] <- ansZ } # 小数部の処理 d <- as.character (round (d, digits)) # 数値誤差があるので丸め d <- rev (as.numeric (strsplit (d, "")[[1]][-(1:2)])) # 0.xxxの0と.は除いておく ansD[1] <- d[1] for (i in seq_len (length (d) - 1)) { ansD <- crossprod (c(.5, 1), d[c (i, i + 1)]) d[i + 1] <- ansD } ans <- ansZ + ansD * .5 cat ("ans = ", ans, "\n") } my2to10 (1101001.1011, digits = 4)
第二章 べき乗と多項式
2進数分解法によるべき乗計算
# 2進数分解法による冪乗計算 myPower <- function (x, n) { binaryN <- strsplit (as.character (my10to2 (n)), "")[[1]] p <- numeric (length (binaryN)) p[1] <- x for (i in seq_len (length (p) - 1)) { p[i + 1] <- p[i] * p[i] } ans <- prod (rev (p) ^ type.convert (binaryN)) return (ans) } myPower (2, 31)
ホーナー法による多項式関数の計算
# 多項式近似の値とその導関数を同時に計算する # ホーナー法 myPolynorminal <- function (x, coef) { poly <- length (coef) y <- coef [1] yprime <- 0 for (i in seq_len (poly - 1)) { yprime <- yprime * x + y y <- y * x + coef [i + 1] } return (list (yp = yprime , y = y)) } coef <- c(-.0064535, .0360885, -.0953294, .1676541, -.2407338, .331799, -.4998741, .9999964, 0) myPolynorminal (1, coef) # 結果の比較 # 計算できてるかはlog (1 + x)の多項式近似でみる x <- seq(0, 2, by = .01) y <- unlist (sapply (x, myPolynorminal, coef)[2, ]) yprime <- unlist (sapply (x, myPolynorminal, coef)[1, ]) ty <- log (1 + x) typ <- 1 / (1 + x) d <- data.frame (y = y, yp = yprime , ty = ty, tpy = typ) matplot(x, d, type = "l", lty = c(1, 5), lwd = 2)
第三章 方程式の解
実数のルートをニュートン法で計算
# *********************************************************************************************** # 実数aのn乗根を求めるニュートン法 myNewton.root.n <- function (a, n, epsilon = 1e-7) { if (a == 0) { cat (n, "roots of ", 0, " = ", 0, "\n") return (0) } if (a == 1) { cat (n, "roots of ", 1, " = ", 1, "\n") return (0) } x.old <- 1 gain <- 1 while (gain > epsilon) { t <- (x.old ^ n - a) / (n * x.old ^ (n - 1)) x.new <- x.old - t x.old <- x.new gain <- abs (t) } return (x.new) } myNewton.root.n (a = 2, n = 2)
実数の逆数をニュートン法で計算
# *********************************************************************************************** # 実数aの逆数を求めるニュートン法 myNewton.reciprocal <- function (a, epsilon = 1e-7) { if (a == 0) { stop ("0はアカンて！！") } x.old <- .1 #初期値はこれくらいにしとかないとすぐ破綻する gain <- 1 while (gain > epsilon) { x.new <- x.old * (2 - a * x.old) gain <- abs (x.new - x.old) cat (x.old, x.new, gain , epsilon , gain > epsilon , "\n") x.old <- x.new } return (x.new) } myNewton.reciprocal (a = 12)
ホーナー法とニュートン法による多項式関数の実数解
# *********************************************************************************************** # 三次方程式をニュートン法で解く # まずホーナー法による多項式の計算をする関数 myPolynorminal3 <- function (x, coef) { poly <- length (coef) y <- coef [1] yprime <- 0 coef.rtn <- numeric (poly); coef.rtn [1] <- y for (i in seq_len (poly - 1)) { yprime <- yprime * x + y y <- y * x + coef [i + 1] coef.rtn[i + 1] <- y } return (list (yp = yprime , y = y, c = coef.rtn)) } # 動作確認 coef <- c(-.0064535, .0360885, -.0953294, .1676541, -.2407338, .331799, -.4998741, .9999964, 0) myPolynorminal3 (1, coef) # ここからが本番 # myNewton.Polynominal <- function (coef, x = NULL, epsilon = 1e-7) { # 多項式の解をホーナー法を使って求める # 初期値を指定するときはxをつかう # エラーがでるときは初期値がおかしいか, 複素数が必要なとき if (is.null(x)) { x.olds <- runif (length (coef) -1) } else { x.olds <- x } ans <- numeric (length (coef) -1) # 解 for (i in seq_along (ans)) { x.old <- x.olds[i] gain <- 1 while (gain > epsilon ) { Y <- myPolynorminal3 (x.old, coef) t <- Y$y / Y$yp x.new <- x.old - t gain <- abs (x.new - x.old) x.old <- x.new } ans [i] <- x.new coef <- Y$c coef <- coef[-length (coef)] # 余りは消しておく # 解じゃなくて, 任意の(x-c)のときはこれをみれば余りがわかる } return (ans) } f <- function (x) 256 * x^3 -130 * x^2 + x + 2 curve (f, -.2, .5, lwd = 2, col = "orange") coef <- c(256, -130, 1, 2) abline (h = 0) abline (v = myNewton.Polynominal (coef = coef)) # 成功している！！！！！！！！！ # ***********************************************************************************************
第四章 連立一次方程式
掃き出し法
# *********************************************************************************************** # 掃き出し法 mySweepingOut <- function (A, b) { nr <- nrow (A) nc <- ncol (A) Ab <- cbind (A, b) # 掃き出し法の部分 index <- 1:nr for (i in seq_len (nr)) { # ピポットが0となっているときの処理 if (Ab[i, i] == 0) { first <- which ((Ab[-(1:i), i] != 0))[1] if (is.na (first)) { Ab[i:nr, 1:nc] <- 0 return (Ab) } else { first <- first + i # 修正2013/11/1 tmp <- Ab[i, ] Ab [i, ] <- Ab[first, ] Ab [first, ] <- tmp } } # ここは普通に計算 Ab[i, ] <- Ab[i, ] / Ab[i, i] for (j in index[-i]) { Ab[j, ] <- Ab[j, ] - Ab[j, i] * Ab[i, ] } } return (Ab) } # 動作確認1 A <- matrix (c (-1, 2, 2, 2, 1, 4, 1, 0, 3), nrow = 3, ncol = 3, byrow = TRUE) b <- c(-3, 7, 7) (mySweepingOut (A, b)) # 動作確2；0となるピポットがあるとき A <- matrix (c (-3, 2, 1, 3, 2, 4, 2, 3, 3, 1, 0, 2, -1, 2, 1, 1), 4, 4, byrow = TRUE) b <- c (1, 0, 2, 3) (mySweepingOut (A, b)) (solve (A, b)) # 結果が同じであることの確認
掃き出し法を使った逆行列
# *********************************************************************************************** # 掃き出し法を用いた逆行列 myInverseMatrix <- function (A) { nc <- ncol (A) nr <- nrow (A) E <- diag (nr) EA <- mySweepingOut (A, E) cE <- EA[, 1:nc] if (all (diag (cE) == 1)) { return (EA [,-(1:nc)]) } else { warning ("出直して来い！") } } A <- matrix (c (1, 2, -1, -1, 1, 2, 1, -1, 1), 3, 3, byrow = TRUE) myInverseMatrix (A) solve (A)
掃き出し法を使った行列式
# *********************************************************************************************** # 掃き出し法を用いた行列式の計算 myDeterminant <- function (A) { nc <- ncol (A) nr <- nrow (A) n <- 0 # 行置換の回数 # 掃き出し法の部分 index <- 1:nr determinant.element <- numeric (nr) for (i in seq_len (nr)) { # ピポットが0となっているときの処理 if (A[i, i] == 0) { first <- which ((Ab[-(1:i), i] != 0))[1] if (is.na (first)) { A[i:nr, 1:nc] <- 0 return (list (Ab, det = 0)) } else { first <- first + i # 修正2013/11/1 n <- n + 1 tmp <- A[i, ] A [i, ] <- A[first, ] A [first, ] <- tmp } } # ここは普通に計算 determinant.element [i] <- A[i, i] A[i, ] <- A[i, ] / A[i, i] for (j in index[-i]) { A[j, ] <- A[j, ] - A[j, i] * A[i, ] } } return (list (A, det = prod (determinant.element) * (- 1) ^ n)) } A <- matrix (c (1, 2, -1, -1, 1, 2, 1, -1, 1), 3, 3, byrow = TRUE) myDeterminant (A) det (A)
ガウスの消去法
# *********************************************************************************************** # ガウスの消去法 myGaussSweep <- function (A, b) { nr <- nrow (A) nc <- ncol (A) Ab <- cbind (A, b) # 消去法の部分 index <- 1:nr for (i in seq_len (nr)) { # ピポットが0となっているときの処理 if (Ab[i, i] == 0) { first <- which ((Ab[-(1:i), i] != 0))[1] if (is.na (first)) { Ab[i:nr, 1:nc] <- 0 return (Ab) } else { first <- first + i # 修正2013/11/1 # 違う行で tmp <- Ab[i, ] Ab [i, ] <- Ab[first, ] Ab [first, ] <- tmp } } # ここは普通に計算 Ab[i, ] <- Ab[i, ] / Ab[i, i] for (j in index[-(1:i)]) { Ab[j, ] <- Ab[j, ] - Ab[j, i] * Ab[i, ] } } return (Ab) } A <- matrix (c (-3, 2, 1, 3, 2, 4, 2, 3, 3, 1, 0, 2, -1, 2, 1, 1), 4, 4, byrow = TRUE) b <- c (1, 0, 2, 3) (myGaussSweep (A, b))
後退代入
# *********************************************************************************************** # ガウスの消去法を使ってからの # 後退代入 # Aは上三角行列 myBackSolve <- function (Ab) { nr <- nrow (Ab) nc <- ncol (Ab) ; ncc <- nc - 1 ans <- numeric (nr) ans [nr] <- Ab[nr, nc] / Ab [nr, nr] for (i in (nr-1):1) { ans [i] <- (Ab[i, nc] - sum (Ab[i, (i + 1):ncc, drop = TRUE] * ans [(i+1):nr])) / Ab[i, i] } return (ans) } myBackSolve (myGaussSweep (A, b)) solve (A, b)
前進代入
# ********************************************************************************************* # 前進代入 myForwardSolve <- function (Ab) { nr <- nrow (Ab) nc <- ncol (Ab) ; ncc <- nc - 1 ans <- numeric (nr) ans [1] <- Ab[1, nc] / Ab [1, 1] for (i in 2:nr) { ans [i] <- (Ab[i, nc] - sum (Ab[i, 1:(i-1), drop = TRUE] * ans [1:(i-1)])) / Ab[i, i] } return (ans) }
ピポットの部分選択(全体的に修正 2013/11/8)
# ********************************************************************************************* # ピポットの部分選択 # ピポットの絶対値が大きいほうから計算したほうがいいらしい myPartPivotSelect <- function (A) { nr <- nrow (A) P <- diag (nr) for (i in 1:(nr-1)) { iP <- diag (nr) max.row.index <- which.max (abs (A[-(1:i), i])) + i if (abs (A[i, i]) < abs (A[max.row.index, i])) { tmp <- A[i, ] A[i, ] <- A[max.row.index, ] A[max.row.index, ] <- tmp tmp <- iP[i, ] iP[i, ] <- iP[max.row.index, ] iP[max.row.index,] <- tmp } P <- crossprod (t (iP), P) } # Aが行置換後の行列 # Pが行置換行列 return (list (mat = A, P = P)) } # 動作確認 A <- matrix (c (1, 2, 3, 2, 4, 8, 2, 3, 8), 3, 3, byrow = TRUE) print (myPartPivotSelect (A)) print (myPartPivotSelect (A)$P %*% A)
ピポットの完全選択
# ********************************************************************************************* # ピポットの完全選択 # 行だけでなく列も含めて考える myAbsolutePivotSelect <- function (A) { nr <- nrow (A) absA <- abs (A) col.index <- 1:ncol (A) row.index <- 1:nrow (A) for (i in 1:(nr-1)) { partA <- abs (A [-(1:i), -(1:i), drop = FALSE]) index <- which.max (partA) nrr <- nrow (partA) cl <- index %/% nrr # 商 rw <- index %% nrr # あまり if (cl == 0) { cl <- cl + 1 } if (rw == 0) { rw <- nrr } if (abs (A[i, i]) < abs (partA[rw, cl])) { # 行の対処 tmp1 <- A[i, ] A[i, ] <- A[i+rw, ] A[i+rw, ] <- tmp1 # 行交換の記録 row.index[i] <- i+rw row.index[i+rw] <- i # 列の対処 tmp2 <- A[, i] A[, i] <- A[, i+cl] A[, i+cl] <- tmp2 # 列交換の記録 col.index[i] <- i+cl col.index[i+cl] <- i } } return (list (mat = A, c.index = col.index, r.index = row.index)) } myAbsolutePivotSelect (A)
LU分解
# ********************************************************************************************* # LU分解クラウト # LDUにするときはLDU ＝TRUE myLUdecomposition <- function (A, LDU = FALSE) { if (!exists ("myBackSolve")) { source ("myBackSolve.r") } if (!exists ("myForwardSolve")) { source ("myFOrwardSolve.r") } nr <- nrow (A) nc <- ncol (A) L <- diag (nr) # ガウスの消去法 index <- 1:nr for (i in seq_len (nr)) { # ピポットが0となっているときの処理 if (A[i, i] == 0) { first <- which ((Ab[-(1:i), i] != 0))[1] if (is.na (first)) { warning ("出直して来い！") } else { first <- first + i # 修正2013/11/1 # tmp1 <- A[i, ] A [i, ] <- A[first, ] A [first, ] <- tmp1 # tmp2 <- L[i, ] L[i, ] <- L[first, ] L[first, ] <- tmp2 } } # 普通に計算 # Lは先に参照するのだ L[i:nr, i] <- A[i:nr, i] A[i, ] <- A[i, ] / A[i, i] for (j in index[-(1:i)]) { A[j, ] <- A[j, ] - A[j, i] * A[i, ] } } # デカイときはこれを使ってもいいかも # library (Matrix) # A <- as (A, "sparseMatrix") # L <- as (L, "sparseMatrix") dimnames (L) <- dimnames (A) if (LDU) { diagonal <- diag (L) L <- sweep (L, 2, c(diagonal), "/") D <- diag (diagonal) dimnames(D) <- dimnames (A) return (list (U = A, L = L, D = D)) } return (list (U = A, L = L)) } A <- matrix (c (-3, 2, 1, 3, 2, 4, 2, 3, 3, 1, 0, 2, -1, 2, 1, 1), 4, 4, byrow = TRUE) rtn <- myLUdecomposition (A, LDU = TRUE) rtn$L %*% rtn$D %*% rtn$U #libraryを使って確かめ library(Matrix) elu <- expand(lu(A)) elu$L %*% elu$U # LU分解を使いつつ連立方程式を解く A <- matrix (c (-3, 2, 1, 3, 2, 4, 2, 3, 3, 1, 0, 2, -1, 2, 1, 1), 4, 4, byrow = TRUE) b <- c (1, 0, 2, 3) LU <- myLUdecomposition (A) # Lの対角成分の積をとれば行列式になる y <- myForwardSolve (cbind (LU$L, b)) x <- myBackSolve (cbind (LU$U, y)) x solve (A, b)
第五章 反復法と収束
データとメモ
A <- matrix (c (4, 1, 0, -2, 1, 4, -1, 1, 0, -1, 4, -1, -2, 1, -1, 4), 4, 4, byrow = TRUE) b <- c(5, -5, 6, -8) solve (A, b) # [1] 1 -1 1 -1 # ********************************************************************************************* # 反復系はちゃんと対角優位にピポットを変えられないと動かいないよ
ヤコビ法
# ********************************************************************************************* # ヤコビ法による連立方程式の解法 # ネットワークや構造力学などの帯行列:疎行列がでてくるときに # ガウスの消去法などよりも有効 # ヤコビ法が収束する条件は特性根が1より小さいこと # -1/log10(|lambda_max|)解ごとに一桁改善する # 実際には固有値はめんどくさいので対角有意かどうかを調べる myJacobiIteration <- function (A, b, epsilon = 1e-15) { nr <- nrow (A) d <- b / diag (A) B <- - sweep (A, 1, diag (A), "/"); diag (B) <- 0 x.old <- runif (nr) gain <- 1 while (gain > epsilon) { x.new <- d + B %*% x.old gain <- max (abs (x.new - x.old)) x.old <- x.new } return (x.new) } myJacobiIteration (A, b)
ガウスサイデル（コツコツ）
# *********************************************************** # ガウスザイデル法その１ myGaussSeidelIteration1 <- function (A, b, epsilon = 1e-15) { nr <- nrow (A) x <- runif (nr) d <- diag (A) gain <- 1 delta <- numeric (nr) while (gain > epsilon ) { for (i in 1:nr) { x.new <- (b[i] - crossprod (A[i, -i], x[-i])) / d[i] delta [i] <- abs (x.new - x[i]) x[i] <- x.new } gain <- max (delta) } return (x) } myGaussSeidelIteration1 (A, b)
ガウスサイデル（行列使いつつ）
# *********************************************************** # ガウスザイデル法その２ myGaussSeidelIteration2 <- function (A, b, epsilon = 1e-15) { nr <- nrow (A) d <- b / diag (A) B <- sweep (A, 1, diag (A), "/"); diag (B) <- 1 C <- B; C[upper.tri (B)] <- 0 D <- -B; D[lower.tri (B, diag = TRUE)] <- 0 x.old <- numeric (nr) gain <- 1 while (gain > epsilon) { y <- d + D %*% x.old x.new <- myForwardSolve (cbind (C, y)) gain <- max (abs (x.new - x.old)) x.old <- x.new } return (x.new) } myGaussSeidelIteration2 (A, b)
第六章 数値積分
メモ
n分割台形則
# ********************************************************************************************* # log 2 を積分で得る # n区間台形則 myTrapezoidalIntegration <- function (FUN, n, lower = 0 , upper = 1) { # nは区間 # n区間なので点はn+1必要 x <- seq (lower, upper, length.out = n + 1) h <- (upper - lower) / n y <- sapply (x, FUN) coef <- c (.5, rep (1, n - 1), .5) Y <- h * crossprod (coef, y) return (Y) } f <- function (x) {1 / (1 + x)} myTrapezoidalIntegration (FUN = f, n = 2)
精度付き台形則
# ********************************************************************************************* # 精度付き台形則 myTrapesoidalIntegration.with.accuracy <- function (FUN, lower = 0, upper = 1, epsilon = 1e-10) { n <- 1 h <- upper - lower I <- h * .5 * (FUN (lower) + FUN (upper)) gain <- 1 while (gain > epsilon) { # 現在の点を求めて, そこからh/2だけズラす. ちょっと非効率かも？ #dx <- seq (lower, upper, length.out = n + 1)[-(n + 1)] + h k <- 1:n dx <- lower + (k - .5) * h D <- .5 * h * sum (sapply (dx, FUN)) gain <- abs (D - I/2) I <- I / 2 + D n <- n * 2 h <- h / 2 } return (I) } f <- function (x) {1 / (1 + x)} myTrapesoidalIntegration.with.accuracy (FUN = f, upper = 1)
2k分割シンプソン則
# *********************************************************************************************** # n区間シンプソン則 mySimpsonIntegration <- function (FUN, n, lower = 0 , upper = 1) { # n区間 if (n%%2 != 0) { warning ("出直して来いb") } # n区間なので点はn+1必要 x <- seq (lower, upper, length.out = n + 1) h <- (upper - lower) / n y <- sapply (x, FUN) coef <- c (1, rep (c (4, 2), length = n - 1), 1) Y <- h * crossprod (coef, y) / 3 return (Y) } f <- function (x) {1 / (1 + x)} mySimpsonIntegration (FUN = f, n = 30)
精度付きシンプソン則
# ********************************************************************************************* # 精度付きシンプソン則 mySimpsonIntegration.with.accuracy <- function (FUN, lower = 0 , upper = 1, epsilon = 1e-7) { n <- 2 h <- (upper - lower) / n f0 <- FUN (lower) f1 <- FUN (h) f2 <- FUN (upper) I <- h * (f0 + 4 * f1 + f2) / 3 J <- 4 * h * f1 / 3 gain <- 1 while (gain > epsilon) { k <- 1:n dx <- lower + (k - .5) * h # Dはひとつ前のJの-1/4倍 D <- - J / 4 J <- 2 * h * sum (sapply (dx, FUN)) / 3 D <- D + J gain <- abs (D - I/2) I <- I / 2 + D n <- n * 2 h <- h / 2 } return (I) } f <- function (x) {1 / (1 + x)} mySimpsonIntegration.with.accuracy (FUN = f, upper = 1, epsilon = 1e-6)
n分割ロンバーグ積分
# *********************************************************************************** # 加速 I(n)にだけでなく, I(n/2)なども組み合わせるとより精度がよい結果を得られる # 台形則を加速させる # ロンバーグ積分 myRombergIntegration <- function (FUN, lower = 0, upper = 1, k) { K <- 1:k N <- 2 ^ (K-1) # 結果を保存するマトリックスをアロケート result.mat <- matrix (0, nrow = k, ncol = k) # 一列目を埋める for (i in K) { result.mat [i, 1] <- myTrapezoidalIntegration (FUN = FUN,lower = lower, upper = upper, n = N[i]) } # 後は順に埋める for (j in 2:k) { # 行 for (u in 2:j) { # 列 p4 <- 4 ** u result.mat[j, u] <- (p4 * result.mat [j, u-1] - result.mat [j-1, u-1]) / (p4 - 1) } } I <- result.mat [k, k] return (I) } f <- function (x) {1 / (1 + x)} myRombergIntegration (FUN = f, k = 20)
精度付きロンバーグ積分
# *********************************************************************************** # 精度付きロンバーグ積分 myRombergIntegration.with.accuracy <- function (FUN, lower = 0, upper = 1, epsilon = 1e-5) { n <- 1 # 分割数 k <- 1 # ステップ数 gain <- 1 # 結果を保存するマトリックスをアロケート result.mat <- matrix (0, nrow = 30, ncol = 30) result.mat [1, 1] <- myTrapezoidalIntegration (FUN = FUN,lower = lower, upper = upper, n = n) I <- result.mat [1, 1] while (gain > epsilon) { k <- k + 1 n <- n * 2 result.mat [k, 1] <- myTrapezoidalIntegration (FUN = FUN, lower = lower, upper = upper, n = n) for (i in 2:k) { p4 <- 4 ** i result.mat[k, i] <- (p4 * result.mat [k, i-1] - result.mat [k-1, i-1]) / (p4 - 1) } I <- result.mat [k, k] gain <- abs (result.mat[k, k] - result.mat[k-1, k-1]) } return (I) } f <- function (x) {1 / (1 + x)} myRombergIntegration.with.accuracy (FUN = f, epsilon = 1e-5)
*1:Dirk Eddelbuettel and Romain Francois (2011). Rcpp: Seamless R and C++ Integration. Journal of Statistical Software, 40(8), 1-18. URL http://www.jstatsoft.org/v40/i08/.
*2:ans W Borchers (2013). pracma: Practical Numerical Math Functions. R package version 1.5.5. http://CRAN.R-project.org/package=pracma
