読者です 読者をやめる 読者になる 読者になる

Rで普通の数値解析その1

塩には必要だが。

あっ, どーも僕です。

来年から

おかげさまで、来年からの食い扶持が決まりました。
特に@beroberoさんは見ず知らずの私に声をかけてくださりありがとうございました。励みになりました。

内定先ではデータ解析とかとはあまり関係ない仕事かもですが、分野としては志望通りで一安心です。

Rで普通の数値解析

近頃というか、前々からのですが、これからのためには普通の数値解析をちゃんと実装したことほうがいいのでわと感じていました。

そこで、今回はRを使ってですが普通の数値解析を勉強しました。今週を使ってコツコツやりました。ただ、ここで実装したのはすべてRで提供されていますのであくまで自分の訓練用です。

それなら、せめてRcpp*1使えって話ですがRcppの教科書は一回しか読んでなくてまだ使いこなせていのでやめました。

教科書はこちら。

数値で学ぶ計算と解析

数値で学ぶ計算と解析

初読は一年くらい前ですが、そのときはただ読んだだけでした。内容のレベルはまえがきによると情報系の学部生向けです。


ちなみに数値解析をRでやるならpracma*2がオススメです。



ってことで、以下には各章ごとに実装したRのプログラムをいくつか載せます。気になる項目があったら見てみてください。一応、全てちゃんと動作しています。ただ、ピポットの選択は自信ないです。(教科書に例題がなかった)また、実装したの次のものです。ところどころ、章をまたいで依存してるかもです。

あと、誤差のことは説明するの大変なので、本買ってください。


第一章 数値と誤差
  • 十進法を二進法に変換
  • 二進法を十進法に変換
第二章 べき乗と多項式
  • 2進数分解法によるべき乗計算
  • ホーナー法による多項式関数の計算
第三章 方程式の解
第四章 連立一次方程式
  • 掃き出し法
  • 掃き出し法を使った逆行列
  • 掃き出し法を使った行列式
  • ガウスの消去法
  • 後退代入
  • 前進代入
  • ピポットの部分選択
  • ピポットの完全選択
  • LU分解
第五章 反復法と収束
  • ヤコビ法
  • ガウスサイデル(コツコツ)
  • ガウスサイデル(行列使いつつ)
第六章 数値積分
  • n分割台形則
  • 精度付き台形則
  • 2k分割シンプソン則
  • 精度付きシンプソン則
  • n分割ロンバーグ積分
  • 精度付きロンバーグ積分


七章は特にやることなかったです。


追記

スクリプトの最新版はこちらから。
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)

ガウスサイデル(コツコツ)


# ***********************************************************
# ガウスザイデル法その1
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)

ガウスサイデル(行列使いつつ)

# ***********************************************************
# ガウスザイデル法その2
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)

第六章 数値積分

メモ

  • dnormを積分してみたらシンプソン則の方が時間かかる
  • シンプソン則は両端の一階微分が0に近いときや, 同じとき収束良くない
  • 次のステップとしては二重指数型積分公式を使う
  • ガウスの積分公式は多項式の厳密な積分値を与える
  • 少ない分割点しか得られないときに有効
  • (つまり, ラグランジュ係数は自力で解く!!)

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