前報の「ロンバーグ積分」には，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
