電脳ヨーグルト

プログラミングの備忘録と私生活を

数値解析の概要と非線形方程式の解を求める二分法について

f:id:at25250410:20180926134237p:plain

今回は数値解析の概要についてざっくり書いてみます。

数値解析とはなにか

数値解析(数値計算)は、方程式の解や関数の値を有限回の演算で「できるだけ早く」「できるだけ正確に」求めるアルゴリズムを模索し,計算結果を検証する学問です。

コンピュータを用いた数値計算では、解析的に解けないもとの問題を
「有限個・有限桁の小数同士の有限回の四則演算」に置きかえて計算します。

数値解析はなんの役に立つのか?

世の中の物理現象の多くは微分方程式によって記述されますが、これらは解析的に(数式としては)解くことができないことが多いです。

例えば、ばね振り子や太陽の周りを回る惑星の運動などがそうです。
数値計算を用いるメリットは、このような解析的に (式変形などをしても手計算では) 解けない問題に対してもおおよその解が得られるという点と、手計算より大規模な問題を扱うことができるという点があります。

二分法を用いて非線形方程式の解を求める

非線形というのは簡単に言うと直線ではないということです。なので直線である一次関数は線形ですが、曲線である二次関数や三次関数は非線形です。

非線形方程式(1次ではない方程式)を解くとは、線形性を持たない(主に曲線) f(x) に対し
f(x) = 0 の等号を成立させる x の値を求めることです。

f:id:at25250410:20180926122718p:plain

上の図のy = f(x)の解は、グラフy = f(x)とx軸が交わる点のx座標(赤い点部分)です。

仮にこのf(x) = X^2 - 2とします。

f(x) = 0を解くと√2という解を出せます。ですが√2は無理数です。
1.41421356......と永遠に続きます。

コンピュータにとってこのような無理数は迷惑なので、二分法を用いて大まかな値を求めます。

以下に二分法を用いた解の算出のアルゴリズムを掲載しています。ちなみにεイプシロンと読みます。
イプシロンは「1より大きい最小の数」と1との差のことで、今回のアルゴリズムで用いるイプシロンは0.00000001のような非常に小さい値だという認識で大丈夫だと思います。

アルゴリズム (二分法)
Step-1 適当な a, b(> a) を選ぶ.
Step-2 もし f(a)f(b) > 0 ならば Step1 に戻る.
Step-3 もし |f(a)| < ε ならば a を解として終了.また,もし

f(b) < ε ならば b を解として終了.

Step-4 c = a+b
2 とし,もし |f(c)| < ε ならば c を解として終了.
Step-5 f(a)f(c) > 0 ならば a = c とし,そうでなければ b = c として
Step4 に戻る.

f:id:at25250410:20180926122718p:plain

解である赤い点を挟むようにして、テキトーにaとbの二点を選択します。次にaとbの中間点cを求めます。この時、f(c)の絶対値がイプシロンを下回った場合はcが解となります。下回らなかった場合はc = aやc = bとなって再び仕切り直しで中間点を求めて新たにcを求めていきます。

このような動作を何十回も繰り返すとcの値は、解に限りなく近い値になります。
f:id:at25250410:20180926133736p:plain

f:id:at25250410:20180926133705p:plain

f(x) = x ^ 2 - 2という方程式の解を二分法によって求めると、上の画像のように繰り返すうちに徐々に正確な値に近づいていることがわかります。

以下のC言語のソースコードを実行すると、二分法を試せますので良かったら参考にしてみてください。

#include <stdio.h>
#include <math.h>
#define EPS   1.0e-10
#define K_MAX     100

void input(double *pa, double *pb);
double f(double x);
double bisection(double a, double b);
void show_each_step(double a, double b, double c, int i);

int main(void) {
    double a, b, ans;

    input(&a,&b);
    ans = bisection(a,b);
    printf("f(%.16f) = %.16f\n",ans,f(ans));

    return 0;
}

void input(double *pa, double *pb) {
    double tmp;

    do {
        printf("input a: ");
        scanf("%lf",pa);
        printf("input b: ");
        scanf("%lf",pb);
    } while (f(*pa) * f(*pb) > 0.0);
    if (*pa > *pb) {
        tmp = *pa;
        *pa = *pb;
        *pb = tmp;
    }
}

double f(double x) {
    return (x * x  - 2);
}

double bisection(double a, double b) {
    double c;
    int i;

    for (i = 1; i <= K_MAX; i++) {
        c = (a + b) / 2.0;
				if (f(c) < EPS && f(c) > 0) break;
        else if (f(a)*f(c) > 0) a = c;
        else if (f(b)*f(c) > 0) b = c;
        show_each_step(a,b,c,i);
    }

    return c;
}

int dgtof(int m) {
    int ans = 1;

    while (m >= 10) {
        ans++;
        m = m / 10;
    }

    return ans;
}

void show_each_step(double a, double b, double c, int i) {
    int dgt = dgtof(K_MAX);

    printf("<%*d> [ %.16f, %.16f ], ",dgt,i,a,b);
    printf("f(%.16f) = % .16f\n",(a+b)/2.0,f(c));
}