数値解析の概要と非線形方程式の解を求める二分法について【エンジニアメモ】

数値解析とはなにか

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

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

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

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

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

上の図の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 に戻る．

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

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

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));
}

数値解析を学ぶのにおすすめの本