【数値情報処理a】第2回代数方程式

前回の課題

【数値情報処理a】第1回の答え

数値計算で生じる誤差

丸め誤差
打ち切り誤差
桁落ち
情報落ち

丸め誤差

浮動小数点数の仮数部の上限の桁数を超えたときに丸めることによって生じる誤差である。

打ち切り誤差

無限に繰り返し行って結果が得られる計算処理を途中で打ち切ることにより生じる誤差である。

例えば、等比数列の和 \[ S_{n}=\sum_{i=1}^{n}(0.5)^{i-1} \] の $n\rightarrow\infty$ の極限を考えよう。

等比数列の和の公式から \[ \lim_{n\rightarrow\infty}S_{n}=\lim_{n\rightarrow\infty}\frac{1-(0.5)^{n}}{1-0.5}=\frac{1}{0.5}=2 \] となることは分かる。

R で $S_{n}$ を計算したらどうなるか見てみよう。無限回の和をとることは不可能なので $n=10$ で近似して計算すると次のようになる。

> s <- 0; for (i in 1 : 10) {s <- s + 0.5 ^ (i - 1)}; print(s)
[1] 1.998047

和を 10 回で打ち切っているため、$2$ にならない。この場合の打ち切り誤差は $2-1.998047=0.001953$ である。

$n$ を増やせば打ち切り誤差は小さくなっていくことが確認できる。

桁落ち

値がほぼ等しく丸め誤差を持つ二つの数値の引き算を行った時に生じる誤差のことである。

有名な例は二次方程式の解の計算で起こる桁落ちである。

$ax^{2}+bx+c=0$ の解は、解の公式より \[ x_{\pm}=\frac{-b\pm\sqrt{b^{2}-4ac}}{2a} \] であるが、$b$ に比べて $a,c$ が非常に小さいときに片方の解が桁落ちを起こしやすい。

$a=b=1, c=10^{-8}$ として片方の解 $x_{+}$ を計算すると、正しくは \[ x_{+}= -1.00000001000\cdots\times 10^{-8} \] となるが、実際に R で計算すると

> a <- 1; b <- 1;c <- 1e-8; print((-b + sqrt(b ^ 2 - 4 * a * c))/(2 * a), digits = 22)
[1] -1.000000010575874398455e-08

となり、小数第10位に誤差が生じる。

この桁落ちを防ぐには、工夫してほぼ等しい数同士の引き算を無くせばよい。 $x_{+}$ は次のように変形できる。 \[ x_{+}=\frac{(-b+\sqrt{b^{2}-4ac})(-b-\sqrt{b^{2}-4ac})}{2a(-b-\sqrt{b^{2}-4ac})}=-\frac{2c}{b+\sqrt{b^{2}-4ac}} \] 変形した式だと、ほぼ等しい数同士の引き算がないので桁落ちを防ぐことができる。

この式で計算すると

> a <- 1; b <- 1;c <- 1e-8; print(-2 * c / (b + sqrt(b ^ 2 - 4 * a * c)), digits = 22)
[1] -1.000000010000000354983e-08

となる。

情報落ち

絶対値が大きい数と小さい数を足したり引いたりしたときに小さい数の情報が落ちてしまうことである。

例えば、$10^{16}$ に $10$ を足してみよう。

> print(1e16 + 10, digits = 22)
[1] 10000000000000010

$10000000000000010$ という結果になり、情報が残っているように見える。

一方、$10^{16}$ に $5$ を足すと情報が落ちてしまうので $10000000000000005$ にならない。

> print(1e16 + 5, digits = 22)
[1] 10000000000000004

$10^{16}$ に $5$ を 2 回足したとしても情報が落ちているので $10000000000000010$ にならない。

> print(1e16 + 5 + 5, digits = 22)
[1] 10000000000000008

代数方程式の解法

代数方程式とは \[ 多項式=0 \] の形で表される方程式である。よくあるのは未知数を $x$ として、代数方程式 \[ f(x)=0 \] の解を求めるという問題である。これは $y=f(x)$ と $y=0$ （$x$ 軸）の交点を見つけるという問題と同じである。

代数方程式の解

例えば、代数方程式が $x$ の 2 次方程式 \[ f(x)=ax^{2}+bx+c=0 \] の場合、その解の公式はよく知られたように \[ x_{\pm}=\frac{-b\pm\sqrt{b^{2}-4ac}}{2a} \] である。

4 次方程式までは代数的解法（有限回の四則演算とべき乗演算）で解くことができて解の公式もあるが、5 次以上や三角関数や指数関数になると代数的解法で解くことが難しくなってくる。そのような場合は、数値計算を使って（近似的に）解くことができる。

代数方程式 $f(x)=0$ において、$x$ にある値を代入して $f(x)=0$ を満たせば、その値が解の一つになる。単純な方法として $x$ を少しずつ変えて解を探す方法があるが、時間がかかって効率が悪い。

効率よく代数方程式の解を求める方法としては、ある仮の解（初期値）から出発して反復しながら解を改良していき、ある条件に達したら反復を終了するのが一般的である。

反復を終了する条件は、理想としては $f(x)=0$ なのだが、浮動小数点数の丸め誤差のために必ずしも $f(x)=0$ になるとは限らない。そこで、$|f(x)|<\delta $ ¹⁾となるか、反復が最大回数に達するかのどちらかの条件を満たした場合に反復を終了するようにする。ここで $\delta$ ²⁾ は十分小さい数で、欲しい解の精度に応じて適切な値に設定する。反復の最大回数を設定するのは反復を無駄に繰り返さないためであり、反復の最大回数に達した場合は十分な精度の解が得られていない可能性が高い。

代数方程式を数値的に解くにはいくつかの方法があるが、解への収束の速さや信頼性に違いがある。ここでは代表的な二分法とニュートン法を扱う。

二分法

ある 2 点間に代数方程式 $f(x)=0$ の解が一つあることが分かっている場合、2 点の中点を使って逐次的に解を求める方法を二分法という。

二分法では以下の手順で解を求める。

2 点 $x_{1}, x_{2}$ を適当に選ぶ。ただし、$f(x_{1})$ と $f(x_{2})$ の符号が異なるように選ぶ。
2 点の中点 $x_{\mathrm{m}}=(x_{1}+x_{2})/2$ を計算する。
$|f(x_{\mathrm{m}})|<\delta$ を満たせば終了する。
$f(x_{\mathrm{m}})$ の符号が $f(x_{1})$ の符号と等しければ $x_{1}=x_{\mathrm{m}}$ とし、そうでなければ $x_{2}=x_{\mathrm{m}}$ とする。
2 に戻る。

はじめの $x_{1}$ と $x_{2}$ は欲しい解が得られるように適切な値に設定しなくてはならない。

$x_{1}$ と $x_{2}$ の間隔は反復ごとに半分になる。 $(1/2)^{10}=1/1024\simeq 10^{-3}$ なので 10 回の反復で精度は約 3 桁良くなる。

二分法は、他の方法に比べて収束の速度は遅いが、解が存在すれば必ず解が得られるというメリットがある。

二分法

二分法の例

例1

関数 \[ f(x)=x^{3}-9x \] において $f(x)=0$ の解を考える。

まず、$y=f(x)$ のグラフを描画してみよう。 R で関数は function(引数) {関数の中身} で定義する。

> f <- function(x) {return(x ^ 3 - 9 * x)}

この関数は f(x) の x に値を入れると x ^ 3 - 9 * x を計算して返す（return）という定義である。

グラフは plot(関数名, 下限, 上限) を用いて描画できる。 abline(a, b) は $y=ax+b$ の直線を描画する関数である。

> plot(f, -5, 5) # y=f(x)のグラフを描画
> abline(0, 0) # y=0 の線を描画

$y=f(x)$ と $y=0$ の交点が $f(x)=0$ の解である。グラフから分かるように $f(x)=0$ の解は三つある。

$f(x)$ は次のように因数分解できる。 \[ f(x)=x(x-3)(x+3) \] これより $f(x)=0$ の解は $x=-3,0,3$ である。この解が二分法で求められるか試してみよう。

bisection.R は二分法で $f(x)=0$ の解の一つを求めるプログラムである。反復の最大回数を 100 回、$\delta=10^{-10}$ としている。 abs() は絶対値を求める関数で、break は for のループを途中で抜けるためのコマンドである。

bisection.R


# 二分法
 
# 関数の定義
f <- function(x) {
  return(x ^ 3 - 9 * x)
}
 
# 反復の回数
iter <- 100
 
# δ
delta <- 1e-10
 
# 2点の初期値
x1 <- 1
x2 <- 10
 
f1 <- f(x1)
f2 <- f(x2)
 
# 反復
for (k in 1 : iter) {
  # 中点
  xm <- (x1 + x2) * 0.5
  fm <- f(xm)
  print(sprintf("%d回目 %.20f", k, xm), quote = FALSE)
 
  # |f(x)|<δ を満たせば終了
  if (abs(fm) < delta) {
    break
  }
 
  # fmとf1の符号が同じか
  if (fm * f1 >= 0.0) {
    x1 <- xm
    f1 <- fm
  } else {
    x2 <- xm
    f2 <- fm
  }
}

2点の初期値を $x_{1}=1$, $x_{2}=10$ としているので、これを実行すると解の一つ $x=3$ に収束していく。ただし、$|f(x_{\mathrm{m}})|<\delta$ を満たした時点で終了するので $x=3$ に近い値で終了する。

[1] 1回目 5.50000000000000000000
[1] 2回目 3.25000000000000000000
[1] 3回目 2.12500000000000000000
︙
省略
︙
[1] 36回目 3.00000000010186340660
[1] 37回目 3.00000000003637978807
[1] 38回目 3.00000000000363797881

他の2つの解に収束させるには初期値を変える必要がある。

例2

関数 \[ f(x)=x^{2}-2 \] において $f(x)=0$ の解を考える。

bisection.R の関数（5行目）を書き換えて解を求めてみよう。

関数（5行目）の書き換え

  return(x ^ 2 - 2)

$f(x)$ は次のように因数分解できる。 \[ f(x)=\left(x+\sqrt{2}\right)\left(x-\sqrt{2}\right) \] $f(x)=0$ の解は $x=-\sqrt{2}, \sqrt{2}$ である（$\sqrt{2}$ は無理数）。

二分法などでこの $f(x)=0$ の解を求めるということは、$\sqrt{2}$ の値を求めることと同じである。

小テスト①

Moodle Server（非公式）で第2回の小テスト①を受験しなさい。

ニュートン法

ニュートン法（もしくはニュートン・ラフソン法）は代数方程式の数値的解法として一般的に使われる方法である。アイザック・ニュートンが考案し、その後、ジョゼフ・ラフソンが改良したものが使われている。

代数方程式 $f(x)=0$ の解に近い適当な $x$ の初期値を $x^{(0)}$ とすると、その点における $y=f(x)$ の接線は \[ y=f'(x^{(0)})(x-x^{(0)})+f(x^{(0)}) \] である。ここで $f'(x)$ は $f(x)$ の微分である。この接線と $x$ 軸との交点を $x^{(1)}$ とすると \[ 0=f'(x^{(0)})(x^{(1)}-x^{(0)})+f(x^{(0)}) \] が成り立つ。この式を $x^{(1)}$ について解くと \[ x^{(1)}=x^{(0)}-\frac{f(x^{(0)})}{f'(x^{(0)})} \] となる。この $x^{(1)}$ を新たな $x^{(0)}$ として同じ操作を繰り返していくと、$x^{(0)}$ は解に近づいていく。

ニュートン法の手順をまとめると以下になる。。

初期値 $x^{(0)}$ を用意する。
$x^{(1)}=x^{(0)}-\frac{f(x^{(0)})}{f'(x^{(0)})}$ を計算して新たな $x^{(0)}$ とする。
$|f(x^{(0)})|<\delta$ を満たせば終了する。
2 に戻る。

この方法では微分 $f'(x)$ が必要であり、微分を解析的に求められば一番よいが、できない場合は数値的に求める方法もある。

ニュートン法は解の近くでは収束が非常に速く、反復ごとに精度の桁が2倍で良くなっていく。例えば、1回目に精度が3桁良くなったとすると、2回目には6桁、3回目には12桁と良くなっていく。

ニュートン法では、解が存在しても反復で収束せず解が得られない場合がある。特に、重解の場合は、解の近くで $f'(x^{(0)})\simeq 0$ になるので、$x^{(1)}$ が発散してしまい収束が悪くなる。

ニュートン法

小テスト②

Moodle Server（非公式）で第2回の小テスト②を受験しなさい。

ニュートン法の例

例1

関数 \[ f(x)=x^{3}-9x \] において $f(x)=0$ の解を考える。

$f(x)$ の微分は \[ f'(x)=3x^{2}-9 \] である。

newton.R はニュートン法で $f(x)=0$ の解の一つを求めるプログラムである。反復の最大回数を 100 回、$\delta=10^{-10}$ としている。

newton.R


# ニュートン法
 
# 関数の定義
f <- function(x) {
  return(x ^ 3 - 9 * x)
}
 
# 関数の微分
fp <- function(x) {
  return(3 * x ^ 2 - 9)
}
 
# 反復の回数
iter <- 100
 
# δ
delta <- 1e-10
 
# 初期値
x0 <- 3.6
 
f0 <- f(x0)
 
# 反復
for (k in 1 : iter) {
  x0 <- x0 - f0 / fp(x0)
  print(sprintf("%d回目 %.20f", k, x0), quote = FALSE)
  f0 <- f(x0)
  # |f(x)|<δ を満たせば終了
  if (abs(f0) < delta) {
    break
  }
}

$x^{(0)}$ の初期値を $x^{(0)}=3.6$ にすると解の一つ $x=3$ に収束していく。この例では、二分法より収束が非常に早い。

[1] 1回目 3.12289156626506025916
[1] 2回目 3.00689294983658106375
[1] 3回目 3.00002362970248048768
[1] 4回目 3.00000000027917623768
[1] 5回目 3.00000000000000000000

他の解に収束させるには $x^{(0)}$ の初期値を変える必要がある。

例2

関数 \[ f(x)=x^{3}-9x+28 \] において $f(x)=0$ の解を考える。これは例1の関数に 28 を足したものである。

$f(x)$ の微分は \[ f'(x)=3x^{2}-9 \] で例1と同じである。

newton.R の関数（5行目）を書き換えて解を求めてみよう。

関数（5行目）の書き換え

  return(x ^ 3 - 9 * x + 28)

$f(x)$ は次のように因数分解できる。 \[ f(x)=(x+4)(x^{2}-4x+7) \] $f(x)=0$ の解は $x=-4$ の一つしかない（虚数解を除く）。

初期値が $x^{(0)}=3.6$ の場合、解 $x=-4$ から遠くて非常に収束が悪い。初期値を色々変えてみると収束性が変わることが分かる。

例3

関数 \[ f(x)=x^{2}-4x+4 \] を考える。

微分は \[ f'(x)=2x-4=2(x-2) \] である。

newton.R の関数（5行目）とその微分（10行目）を書き換えて解を求めてみよう。

関数（5行目）の書き換え

  return(x ^ 2 - 4 * x + 4)

微分（10行目）の書き換え

  return(2 * x - 4)

$f(x)$ は次のように因数分解できる。 \[ f(x)=(x-2)^{2} \] $f(x)=0$ の解 $x=2$ は重解で $f^{\prime}(2)=0$ になってしまうので、ニュートン法では収束が悪くなる。

課題

Moodle Server（非公式）で第2回の課題を行いなさい。

締め切り：2026年4月22日（水）20時

¹⁾

$|f(x)|$ は $f(x)$ の絶対値。

²⁾

ギリシャ文字のデルタ

【数値情報処理a】第2回 代数方程式

前回の課題

数値計算で生じる誤差

丸め誤差

打ち切り誤差

桁落ち

情報落ち

代数方程式の解法

二分法

二分法の例

例1

例2

小テスト①

ニュートン法

小テスト②

ニュートン法の例

例1

例2

例3

課題

【数値情報処理a】第2回代数方程式