プログラミングをするときに知っておきたい数値計算のアルゴリズム

  • 31
    いいね
  • 0
    コメント

今回は主に数値解析のプログラミングをするときに知っておくと便利なアルゴリズムを紹介します。なるべくサンプルソースを付けるようにしますが、その分長くなってしまったのはご愛敬1

多項式の計算

Horner 法

理論

例えば、三次式 $f(x) = a_0 + a_1 x + a_2 x^2 + a_3 x^3$ をそのまま計算すると乗算 6 回、加算 3 回が必要だが、

\begin{align}
f(x) &= a_0 + a_1 x + (a_2 + a_3 x)x^2 \\
     &= a_0 + (a_1 + (a_2 + a_3 x)x)x
\end{align}

と変形してから計算すれば、乗算 3 回、加算 3 回となる。多項式の次数が高いほど効果が高い。

サンプルソース

double horner(double a[], double x) {
    double f = 0;
    for(int i = a.length - 1; i >= 0; i--) {
        f = a[i] + f * x;
    }
    return f;
}

非線形方程式を解く

二分法

理論

$f(x)$ は区間 $[a, b]$ で連続かつ $f(a)f(b) < 0$ とする。このとき中間値の定理により $a < c < b$ で $f(c) = 0$ となるものが存在するが、この $c$ を以下のアルゴリズムによって見つけるのが二分法である。簡単のため $f(a) < 0, f(b) > 0$ とする。

  1. $a_1 = a, b_1 = b$ と置く。
  2. $a_{n+1}$ と $b_{n+1}$ を以下のアルゴリズムで決定する。
    • $\frac{a_n + b_n}{2} > 0$ ならば $a_{n + 1} = a_n, b_{n + 1} = \frac{a_n + b_n}{2}.$
    • そうでなければ $a_{n + 1} = \frac{a_n + b_n}{2}, b_{n + 1} = b_n.$

このとき $c = \lim_{n \to \infty} a_n = \lim_{n \to \infty} b_n$ とすれば $f(c) = 0$ である。

サンプルソース

import java.util.function.DoubleUnaryOperator;

public class Equation {
    static double ERROR = 1e-12; /* Error */

    // Bisection method
    double bisection(DoubleUnaryOperator f, double a, double b) {
        double ta = a;
        double tb = b;
        while(tb - ta > ERROR) {
            double tmp = (ta + tb) / 2;
            if(f.applyAsDouble(tmp) * f.applyAsDouble(tb) > 0) {
                tb = tmp;
            } else {
                ta = tmp;
            }
        }
        return (ta + tb) / 2;
    }
}
public static void main(String args[]) {
    DoubleUnaryOperator f = x -> x * x - 2;
    Equation solver = new Equation();
    System.out.println(solver.bisection(f, 1, 2));
}

面倒なので区間が不正でないかのチェックは省略してます。Java8 になって関数をまるっと引数で渡せるのはいいですね。

Newton 法

理論

$f(x)$ は 2 階微分可能で、以下の条件を満たすとする。

  • $f(a) < 0, f(b) > 0.$
  • $f'(x) > 0, f''(x) > 0 (a \leq x \leq b).$

このとき

  • $c_1 = b.$
  • $c_{n + 1} = c_n - \frac{f(c_n)}{f'(c_n)} (n \geq 1)$

とすれば、$c = \lim_{n \to \infty} c_n$ が存在して $f(c) = 0.$
$$|g'(c)| = \left|\frac{f(c)f''(c)}{(f'(c))^2}\right| = 0 < 1$$
のため、非常に収束が早い。

サンプルソース

import java.util.function.DoubleUnaryOperator;

public class Equation {
    static double ERROR = 1e-12; /* Error */

    // Newton method
    double newton(DoubleUnaryOperator f, DoubleUnaryOperator df, double init) {
        double c = init;
        while(Math.abs(f.applyAsDouble(c)) > ERROR) {
            c -= f.applyAsDouble(c) / df.applyAsDouble(c);
        }
        return c;
    }
}
public static void main(String args[]) {
    DoubleUnaryOperator f = x -> x * x - 2;
    DoubleUnaryOperator df = x -> 2 * x;
    Equation solver = new Equation();
    System.out.println(solver.newton(f, df, 2));
}

数値積分

数値積分は、積分区間を分割(計算の都合上等分するのが慣例)して、その区間の積分を何らかの形で近似計算し、合計することによって求める方法。今回は分かりやすく有界区間上の積分に限定する。

以下、$f(x)$ は区間 $[a, b]$ で連続な関数とする。

台形公式

理論

区間 $[a, b]$ を $n$ 等分し、$h = \frac{b - a}{n}$ とする。
$a_i = a + ih (i = 0, 1, \dots, n)$ とする。$a_0 = a, a_n = b$ である。
区間 $[a_i, a_{i + 1}]$ $(i = 0, 1, \dots, n - 1)$ において
$$\int_{a_i}^{a_{i + 1}}f(x) dx$$

$$\frac{h}{2}(f(a_i) + f(a_{i + 1}))$$
で近似することにより
$$\int_a^b f(x)dx$$

$$\frac{h}{2}\sum_{i = 0}^{n = 1} (f(a_i) + f(a_{i + 1}))$$
で近似できる。

サンプルソース

static double trapezoid(DoubleUnaryOperator f, double a, double b, int step) {
    double h = (b - a) / step;
    double s = 0;
    for(int i = 0; i < step; i++) {
        s += (f.applyAsDouble(a + i * h)
                + f.applyAsDouble(a + (i + 1) * h))
                * h / 2;
    }
    return s;
}

Simpson の公式

理論

区間 $[a, b]$ を $2n$ 等分し、$h = \frac{b - a}{2n}$ とする。
$a_i = a + ih (i = 0, 1, \dots, 2n)$ とする。$a_0 = a, a_{2n} = b$ である。
区間 $[a_{2i}, a_{2i + 2}]$ $(i = 0, 1, \dots, n - 1)$ において
$$\int_{a_{2i}}^{a_{2i + 2}}f(x) dx$$
を、$f(x)$ を三点 $(a_{2i}, f(a_{2i})), (a_{2i + 1}, f(a_{2i + 1})), (a_{2i + 2}, f(a_{2i + 2}))$ を通る放物線で近似することにより
$$\frac{h}{3}(f(a_{2i}) + 4f(a_{2i + 1}) + f(a_{2i + 2}))$$
で近似することが出来る。これにより
$$\int_a^b f(x)dx$$

$$\frac{h}{3}\sum_{i = 0}^{n - 1} (f(a_{2i}) + 4f(a_{2i + 1}) + f(a_{2i + 2}))$$
で近似できる。

サンプルソース

import java.util.function.DoubleUnaryOperator;

public class Integral {
    public static double simpson(DoubleUnaryOperator f, double a, double b, int step) {
        double h = (b - a) / step;
        double s = 0;
        for(int i = 0; i < step / 2; i++) {
        s += (f.applyAsDouble(a + 2 * i * h)
                    + 4 * f.applyAsDouble(a + (2 * i + 1) * h)
                    + f.applyAsDouble(a + (2 * i + 2) * h))
                    * h / 3;
        }
        return s;
    }
}

実際に
$$\int_1^2 \frac{1}{x} dx = \log 2$$
を計算させる。前回結果との差が EPS (適宜定義せよ)より小さくなるまで、刻み幅を 10 ずつ増やしていく。

public static void main(String args[]) {
    DoubleUnaryOperator f = x -> 1 / x;
    double a = 1, b = 2;
    double s;
    s = Integral.simpson(f, a, b, 1);
    int n = 10;
    while(true) {
        double tmp = Integral.simpson(f, a, b, n);
        if(Math.abs(s - tmp) < EPS) break;
        s = tmp;
        n += 10;
    }
    System.out.println(s);
}

常微分方程式を解く

一階の常微分方程式の初期値問題
$$\boldsymbol{x}' = \boldsymbol{f}(t, \boldsymbol{x}) (\boldsymbol{x} = \boldsymbol{x}(t)), \boldsymbol{x}(t_0) = \boldsymbol{x}_0$$
を数値的に解析する方法。1 段法と多段法があるが、今回は分かりやすい 1 段法のみ紹介する。

Euler 法

理論

$$\boldsymbol{x}'(t) = \lim_{h \to 0} \frac{\boldsymbol{x}(t + h) - \boldsymbol{x}(t)}{h}$$

$$\boldsymbol{x}'_i = \frac{\boldsymbol{x}_{i + 1} - \boldsymbol{x}_i}{h}$$
で近似する。ただし $\boldsymbol{x}_i = \boldsymbol{x}(t_i), t_i = t + ih.$ ここで元の微分方程式により $\boldsymbol{x}_i' = \boldsymbol{f}(t_i, \boldsymbol{x}_i)$ であるから
$$\boldsymbol{x}_{i + 1} = \boldsymbol{x}_i + h \boldsymbol{f}(t_i, \boldsymbol{x}_i)$$
で逐次解曲線を近似できる。ただし誤差は刻み幅に比例するので、刻み幅をかなり細かくしないといけないうえ、それに伴う計算量の増加により丸め誤差が蓄積しやすい弱点がある。

サンプルソース

なし(精度が悪く、あまり使用されるアルゴリズムでないため)

Runge-Kutta 法

理論

これに関しては、次の資料が詳しい。
「PDF形式の数学ノート : MATHEMATICS.PDF」内の「4次のRunge-Kutta法」

結論のみ述べると

\begin{align}
\boldsymbol{x}_{i + 1} &= \boldsymbol{x}_i + \frac{h}{6} (\boldsymbol{k}_1 + 2\boldsymbol{k}_2 + 2\boldsymbol{k}_3 + \boldsymbol{k}_4), \\
\boldsymbol{k}_1 &= \boldsymbol{f}(t_i, \boldsymbol{x}_i), \\
\boldsymbol{k}_2 &= \boldsymbol{f}(t_i + \frac{h}{2}, \boldsymbol{x}_i + \frac{k_1}{2}), \\
\boldsymbol{k}_3 &= \boldsymbol{f}(t_i + \frac{h}{2}, \boldsymbol{x}_i + \frac{k_2}{2}), \\
\boldsymbol{k}_4 &= \boldsymbol{f}(t_i + h, \boldsymbol{x}_i + k_3).
\end{align}

となる。誤差は刻み幅の 4 乗に比例するため、精度は Euler 法よりもかなり上がる。

サンプルソース

$x' = \sin t + \cos x, x(0) = 0$ を数値的に解く。

import java.io.BufferedWriter;
import java.io.File;
import java.io.FileWriter;
import java.io.IOException;
import java.io.PrintWriter;
import java.util.function.DoubleBinaryOperator;

public class RungeKutta {
    public static void main(String args[]) {
        File file = new File("solve.txt");
        double t = 0;
        double x = 0;
        double h = 0.01;
        DoubleBinaryOperator f = (a, b) -> Math.sin(a) + Math.cos(b);
        PrintWriter pw = null;
        try {
            pw = new PrintWriter(new BufferedWriter(new FileWriter(file)));
            for(int i = 0; i <= 200; i++) {
                pw.printf("%.8f %.8f\n", t, x);
                double k1 = f.applyAsDouble(t, x);
                double k2 = f.applyAsDouble(t + h / 2, x + h * k1 / 2);
                double k3 = f.applyAsDouble(t + h / 2, x + h * k2 / 2);
                double k4 = f.applyAsDouble(t + h, x + h * k3);
                t += h;
                x += h * (k1 + 2 * k2 + 2 * k3 + k4) / 6;
            }
        } catch (IOException e) {
            e.printStackTrace();
        } finally {
            if(pw != null) pw.close();
        }
    }
}

出力ファイルは Gnuplot 等で表示させてみると良い。少し手直しして、CSV 形式で出力して、Excel でグラフを表示させても良い。


  1. サンプルソースを付けたいがために久々に Java 書いたら若干心が荒んだのはここだけの話(苦笑)。 

この投稿は 数学 Advent Calendar 201613日目の記事です。