未分類
Romberg積分; Romberg Integrate

概要

---

Romberg積分とはRichardsonの補外法を用いた加速収束により、積分値を台形公式に比べより正確に求める数値積分法の一つ。
Richardsonの補外法は収束することが保証されている数列の収束値を求めることができる。
区間を細分化したときその積分値が収束する関数の積分値は求めることができるが、積分区間に特異点を含むときや、ワイエルシュトラス関数のようなフラクタル性を含むなど微分不連続な関数の積分値は求めることができない。
\(y = \sqrt{1-x^2} \)を\([0, 1]\)で積分する場合など積分区間に微分不可能な点を含む場合は正確な値を求めることが難しい。

またRichardsonの補外法は曲線弧長の算出などにも応用ができる。



Richardsonの補外法では以下のように\(i\)番目の収束値\(v_{i, i}\)を求める。
\(\quad \displaystyle v_{i, j} = \frac{4^j v_{i, j - 1} - v_{i - 1, j - 1}}{4^j - 1} \quad i, j \geq 1 \)

ソースコード

---

namespace Intergral {

    /// <summary>Romberg積分</summary>
    public static class RombergIntegral {

        /// <summary>積分</summary>
        /// <param name="func">対象の関数</param>
        /// <param name="a">下端</param>
        /// <param name="b">上端</param>
        /// <param name="precision_level">精度レベル</param>
        /// <remarks>有効桁数14-15桁 精度レベルは大きすぎるとかえって誤差が大きくなる</remarks>
        public static double Integrate(Func<double, double> func, double a, double b, int precision_level = 10) {
            if(func == null) {
                throw new ArgumentNullException(nameof(func));
            }
            if(!(a <= b)) {
                throw new ArgumentException($"{nameof(a)},{nameof(b)}");
            }
            if(precision_level < 1 || precision_level > 16) {
                throw new ArgumentException(nameof(precision_level));
            }

            int max_div = 1 << precision_level;
            double h = b - a, min_h = (b - a) / max_div;
            double[] v = new double[max_div + 1];
            RichardsonExtrapolation conv = new RichardsonExtrapolation();

            for(int i = 0; i <= max_div; i++) {
                v[i] = func(a + i * min_h);
            }

            double t = h * (v[0] + v[max_div]) / 2, new_t;
            conv.Inject(t);

            for(int s = max_div; s > 1; s /= 2) {
                new_t = 0;
                for(int i = s / 2; i < max_div; i += s) {
                    new_t += v[i];
                }

                h /= 2;
                t = t / 2 + h * new_t;

                conv.Inject(t);
            }

            return conv.ConvergenceValue;
        }
    }

    /// <summary>Richardsonの補外法</summary>
    public class RichardsonExtrapolation {
        readonly List<double[]> values = new List<double[]>();

        /// <summary>収束数列値を挿入</summary>
        public void Inject(double new_value) {
            if(SeriesCount <= 0) {
                values.Add(new double[] { new_value });
                return;
            }

            double[] t = values[SeriesCount - 1], t_next = new double[SeriesCount + 1];

            t_next[0] = new_value;

            double r = 1;
            for(int i = 1; i <= SeriesCount; i++) {
                r *= 4;
                t_next[i] = t_next[i - 1] + (t_next[i - 1] - t[i - 1]) / (r - 1);
            }

            values.Add(t_next);
        }

        /// <summary>補外数列</summary>
        public IEnumerable<double> Series {
            get {
                for(int i = 0; i < values.Count; i++) {
                    yield return values[i][i];
                }
            }
        }

        /// <summary>収束推定値</summary>
        public double ConvergenceValue {
            get {
                if(SeriesCount <= 0) {
                    throw new InvalidOperationException();
                }

                return values[SeriesCount - 1][SeriesCount - 1];
            }
        }

        /// <summary>補外数列の要素数</summary>
        public int SeriesCount => values.Count;

        /// <summary>推定最大誤差</summary>
        public double Epsilon {
            get {
                if(SeriesCount <= 1) {
                    throw new InvalidOperationException();
                }

                return Math.Abs(values[SeriesCount - 1][SeriesCount - 1] - values[SeriesCount - 2][SeriesCount - 2]);
            }
        }
    }
}

単体テスト

---

namespace Intergral.Tests {
    [TestClass()]
    public class RombergIntegralTests {
        [TestMethod()]
        public void IntegrateTest() {
            Func<double, double> f = (x) => Math.Sqrt(1 - x * x);

            double v = RombergIntegral.Integrate(f, 0, Math.Sqrt(2) / 2);

            Assert.AreEqual(v, (Math.PI + 2) / 8, 1e-15);
        }
    }
}