cms.integ_romberg のソースコード

import numpy as np
from math import exp

"""
  Numerial integration by Romberg extrapolation
"""

# define function to be integrated


def func(x):
    return exp(x)


# define analytical integration of f(x)


def F(x):
    return exp(x)


#===================
# parameters
#===================
nhmax = 15
eps   = 1.0e-6
xmin, xmax = 1.0, 2.0

#===================
# main routine
#===================
print("Numerical integration using by Romberg extrapolation method")

fx = func(xmin)
Fx = F(xmax) - F(xmin)
print("")
print("eps for convergnece: {}".format(eps))
print("Analytical values:")
print("  f({})={}".format(xmin, fx))
print("  integ(f)[{}, {}]={}".format(xmin, xmax, Fx))
print("")

# create 2D array by specifing the size to be nhmax * nhmax
S = []
for i in range(nhmax):
  S.append([])
  for j in range(nhmax):
    S[i].append([])

n = 1
h0 = xmax - xmin
y0 = func(xmin)
y1 = func(xmax)
S[0][0] = (y0 + y1) / 2.0 * h0;
print("S00 = {}".format(S[0][0]))

kfin = 0
for k in range(1, nhmax+1):
    n *= 2
    h = h0 / n

# 台形則の計算
    S[k][0] = 0.0;
    for i in range(n):
        xi  = xmin + i * h;
        xi1 = xi + h;
        yi  = func(xi);
        yi1 = func(xi1);
        S[k][0] += (yi + yi1) / 2.0 * h;
        print("  S[{}][0] = {}".format(k, S[k][0]))

    for d in range(1, k+1):
       S[k][d] = (4.0**d * S[k][d-1] - S[k-1][d-1]) / (4.0**d - 1.0);
       print("    S[{}][{}] = {}".format(k, d, S[k][d]))

    diff = S[k][k] - S[k-1][k-1]
    print("  diff = S[{}][{}] - S[{}][{}] = {} - {} = {}"
          .format(k, k, k-1, k-1, S[k][k], S[k-1][k-1], diff))
    if abs(diff) < eps:
        break


"""
for m in range(1, nhmax+1):
    h0 *= 0.5
    S[0][m] = (func(x0+h0) - func(x0-h0)) / 2.0 / h0
    print("S0{}: {}".format(m, D[0][m]))
    for k in range(m-1, -1, -1):
        m1 = m - k
        S[m1][k] = (4.0**m1 * D[m1-1][k+1] - D[m1-1][k]) / (4.0**m1 - 1)
        print("D{}{}: {}".format(m1, k, D[m1][k]))
        break
"""

"""
	my $kfin;
	for(my $k = 1 ; $k < $nMax ; $k++) {
		$n *= 2;
print "k=$k (n=$n)\n" if($IsPrint);
		$h = ($x1 - $x0) / $n;

# 台形則の計算
		$S[$k][0] = 0.0;
		for(my $i = 0 ; $i < $n ; $i++) {
			my $xi  = $x0 + $i * $h;
#print "i=$i, xi=$xi\n";
			my $xi1 = $xi + $h;
			my $yi  = &$pFunc($xi);
			my $yi1 = &$pFunc($xi1);
			$S[$k][0] += ($yi + $yi1) / 2.0 * $h;
		}
print "  S[$k][0] = $S[$k][0]\n" if($IsPrint);

		my $d4 = 4.0;
		for(my $d = 1 ; $d <= $k ; $d++) {
#print "    d=$d\n";
			$S[$k][$d] = ($d4 * $S[$k][$d-1] - $S[$k-1][$d-1]) / ($d4 - 1.0);
print "    S[$k][$d] = $S[$k][$d]\n" if($IsPrint);

			$d4 *= 4.0;
		}
		
		my $delta = abs($S[$k][$k] - $S[$k-1][$k-1]);
		if($delta < $EPS) {
			$kfin = $k;
			last;
		}
	}

	my $mess;
	if(!defined $kfin) {
		$mess = "Not converged";
		$kfin = $nMax-1;
	}
	my $Sfin = $S[$kfin][$kfin];
	my $diff = ($kfin >= 1)? $Sfin - $S[$kfin-1][$kfin-1] : undef;

	return ($Sfin, $n, $kfin, $diff, $mess);
}

"""