アルキメデスの方法による円周率の近似値

直径 1 の円の外接正 n 角形の周長を Pn，内接正 n 角形の周長を pn とすれば，

n = 3 * 2 ^ k, k = 0, 1, 2, ...

P2n = 2 * Pn * pn

p2n = sqrt(P2n * pn)

ということで，以下のプログラム

import scipy as sp

def Pi_Archimedes():
    print("{0:>2s} {1:>9s} {2:>17s} {3:>17s}".format("k", "n", "Pn", "pn"))
    k = 0 # 正三角形から始める
    Pn = 3 * sp.sqrt(3)
    pn = Pn / 2
    while True:
        n = 3 * 2 ** k
        print("{0:2d} {1:9d} {2:.15f} {3:.15f}".format(k, n, Pn, pn))
        if Pn == pn:
            return Pn
        k += 1
        P2n = 2 * Pn * pn / (Pn + pn)
        p2n = sp.sqrt(P2n * pn)
        Pn, pn = P2n, p2n

Pi_Archimedes()

 k         n                Pn                pn
 0         3 5.196152422706632 2.598076211353316
 1         6 3.464101615137754 3.000000000000000
 2        12 3.215390309173473 3.105828541230249
 3        24 3.159659942097500 3.132628613281238
 4        48 3.146086215131435 3.139350203046867
 5        96 3.142714599645368 3.141031950890509
 6       192 3.141873049979824 3.141452472285462
 7       384 3.141662747056849 3.141557607911857
 8       768 3.141610176604690 3.141583892148318
 9      1536 3.141597034321526 3.141590463228050
10      3072 3.141593748771351 3.141592105999271
11      6144 3.141592927385096 3.141592516692157
12     12288 3.141592722038614 3.141592619365384
13     24576 3.141592670701998 3.141592645033691
14     49152 3.141592657867844 3.141592651450767
15     98304 3.141592654659306 3.141592653055036
16    196608 3.141592653857171 3.141592653456104
17    393216 3.141592653656637 3.141592653556370
18    786432 3.141592653606503 3.141592653581437
19   1572864 3.141592653593970 3.141592653587703
20   3145728 3.141592653590837 3.141592653589270
21   6291456 3.141592653590053 3.141592653589661
22  12582912 3.141592653589857 3.141592653589759
23  25165824 3.141592653589808 3.141592653589784
24  50331648 3.141592653589796 3.141592653589790
25 100663296 3.141592653589793 3.141592653589791
26 201326592 3.141592653589792 3.141592653589792
27 402653184 3.141592653589792 3.141592653589792

3.141592653589792

正 4 億多角形ですと！

 

R によるプログラムもほとんど同じで，以下のように

Pi.Archimedes = function() {
    Pn = 3 * sqrt(3)
    pn = Pn / 2
    while (TRUE) {
        if (Pn == pn) {
            return(Pn)
        }
        tmp = 2 * Pn * pn /  (Pn + pn)
        pn = sqrt(tmp * pn)
        Pn = tmp
    }
}

> options(digits=16)
> Pi.Archimedes()
[1] 3.141592653589792

追記

Julia では

function piarchimedes()
    Pn = 3sqrt(3)
    pn = Pn / 2
    while true
        Pn == pn && return Pn
        tmp = 2Pn * pn /  (Pn + pn)
        pn = sqrt(tmp * pn)
        Pn = tmp
    end
end

piarchimedes() # 3.141592653589792

ロンバーグ積分（その2）

前報の「ロンバーグ積分」には，C で書かれたプログラムを Python に書き換えたものを示した。両言語は共に「0 オリジン」

https://jp.mathworks.com/matlabcentral/fileexchange/34-rombint-m には，MATLAB で書かれたプログラムがある。

# ROMBINT     Numerical evaluation of an integral using the Romberg method.
#
#   Q = rombint(F, A, B) approximates the integral of F(X) from A to B to
#   within a relative error of 1e-10 using Romberg's method of integration.
#   F is a function.  The function must return a vector of output values if
#   a vector of input values is given.
#
#   Q = rombint(F, A, B, DECIMALDIGITS) integrates with accuracy 10^{-DECIMALDIGITS}.
#   ---------------------------------------------------------
#   Author: Martin Kacenak,
#           Department of Applied Informatics and Process Control,
#           Faculty of BERG, Technical University of Kosice,
#           B.Nemcovej 3, 04200 Kosice, Slovak Republic
#   E-mail: ma.kac@post.cz
#   Date:   posted February 2001, updated June 2006.

MATLAB は「1 オリジン」なので，同じく「1 オリジン」である R で書き換えるのが面倒がない。
見かけはちょっと違うが，ほぼ同じ長さのプログラムである。

rombint = function(funfcn, a, b, decdigs = 10) {
  rom = matrix(0, 2, decdigs)
  romall = numeric(2 ^ (decdigs - 1) + 1)
  h = b - a
  romall = funfcn(seq(a, b, by = h / 2 ^ (decdigs - 1)))
  rom[1, 1] = h * (romall[1] + romall[length(romall)]) / 2
  for (i in 2:decdigs) {
    st = 2 ^ (decdigs - i + 1)
    rom[2, 1] = (rom[1, 1] + h * sum(romall[seq(st / 2 + 1, 2 ^ (decdigs - 1), st)])) / 2
    for (k in 1:(i - 1)) {
      rom[2, k + 1] = ((4 ^ k) * rom[2, k] - rom[1, k]) / (4 ^ k - 1)
    }
    rom[1,] = rom[2, ]
    h = h / 2
  }
  return(rom[1, decdigs])
}

Pi = function(x) {
  return(4/(1+x^2))
}

> rombint(Pi, 0, 1)
[1] 3.141592653589793

> pi - rombint(Pi, 0, 1)
[1] 4.440892098500626e-16

ロンバーグ積分

ロンバーグ積分とは，まず合成台形公式で近似し、次に Richardson の補外法で近似をより正確にする。

Romberg's method

に，C で書かれたプログラムがある。

これを，Python に書き換えると以下のようになる。

for をベクトル計算にするなどにより，ずいぶん短くなる？

import scipy as sp

def Romberg(f, a, b, max_steps=1000, acc=1e-14):
    Rp = sp.empty(max_steps)  # previous row
    Rc = sp.empty(max_steps)  # current row
    h = b - a
    Rp[0] = (f(a) + f(b)) * h / 2
    for i in range(1, max_steps):
        h /= 2
        Rc[0] = h * sum(f(sp.arange(1, 2 ** i, 2) * h)) + Rp[0] / 2  # R(i,0)
        for j in range(1, i + 1):
            n_k = 4 ** j
            Rc[j] = (n_k * Rc[j - 1] - Rp[j - 1]) / (n_k - 1)  # R(i,j)
        if i > 1 and abs(Rp[i - 1] - Rc[i]) < acc:
            return Rc[i - 1]
        Rp, Rc = Rc, Rp # swap Rn and Rc
    return Rp[max_steps - 1]  # best guess

def Pi(x):
    return 4 / (1 + x**2)

a = Romberg(Pi, 0, 1, max_steps=9, acc=1e-14)
a　# 3.141592653589794
sp.pi - a  # -8.881784197001252e-16

積分について

積分は，台形公式，シンプソンの公式のほかにブールの公式がある。

目的とする関数は以下のものを例示する

R

Pi = function(x) {
  return(4/(1+x^2))
}

Python

def Pi(x):
    return 4/(1+x**2)

Trapezoidal rule

プログラムは以下のようになる

R

Trapezoidal = function(f, lo, hi, n = 1000) {
  x = seq(lo, hi, length=n+1)
  w = diff(x)[1]
  return(w*(f(x[1])/2 + sum(f(x[2:n])) + f(x[n+1])/2))
}

> pi - Trapezoidal(Pi, 0, 1)
[1] 1.666666662458738e-07

Python

def Trapezoidal(f, lo, hi, n = 1000):
    x = sp.linspace(lo, hi, n+1)
    w = sp.diff(x)[0]
    return w*(f(x[0])/2 + sum(f(x[1:n])) + f(x[n])/2)

sp.pi - Trapezoidal(Pi, 0, 1) # 1.6666666891040904e-07

Simpson's rule

合成シンプソン則のプログラムは以下のようになる

R

Simpson = function(f, lo, hi, n = 1000) {
  x = seq(lo, hi, length=n+1)
  w = diff(x)[1]
  return(w*(f(x[1])+4*sum(f(x[seq(2, n, by = 2)]))+2*sum(f(x[seq(3, n, by = 2)]))+f(x[n+1]))/3)
}

> pi - Simpson(Pi, 0, 1)
[1] 0

Python

def Simpson(f, lo, hi, n = 1000):
    x = sp.linspace(lo, hi, n+1)
    w = sp.diff(x)[0]
    return w*(f(x[0])+4*sum(f(x[1:n:2]))+2*sum(f(x[2:n:2]))+f(x[n]))/3

sp.pi - Simpson(Pi, 0, 1) # -3.1086244689504383e-15

Boole's rule

プログラムは以下のようになる

R

Boole = function (f, lo, hi, n = 1000) {
  n = ceiling(n / 4) * 4
  x = seq(lo, hi, length=n+1)
  w = diff(x)[1]
  return(2*w*(
      7*sum(f(x[seq(1, n, by=4)]))
    +32*sum(f(x[seq(2, n, by=4)]))
    +12*sum(f(x[seq(3, n, by=4)]))
    +32*sum(f(x[seq(4, n, by=4)]))
    + 7*sum(f(x[seq(5, n+1, by=4)])))/45)
}

> pi - Boole(Pi, 0, 1)
[1] 0

Python

def Boole(f, lo, hi, n = 1000):
    n = int(sp.ceil(n / 4) * 4)
    x = sp.linspace(lo, hi, n+1)
    w = sp.diff(x)[0]
    return 2*w*(
        7*sum(f(x[0:n:4]))
        +32*sum(f(x[1:n:4]))
        +12*sum(f(x[2:n:4]))
        +32*sum(f(x[3:n:4]))
        +7*sum(f(x[4:n+1:4])))/45

sp.pi - Boole(Pi, 0, 1) # 0.0

 

ブールの公式はなかなかのものであることがわかった。