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

===== 前回の課題 =====

[[ans2]]

===== 代数方程式の解法 =====

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

例えば、代数方程式が $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$ を少しずつ変えて解を探す方法があるが、時間がかかって効率が悪い。

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

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

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

===== 二分法 =====
ある 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 桁良くなる。

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

{{ bisection-method.png?nolink |二分法}}


===== 二分法の例 =====
==== 例1 ====

関数
\[
f(x)=x^{3}-9x
\]
を考える。

まず、$y=f(x)$ のグラフを描画してみよう。
R でグラフは ''plot(関数名, 下限, 上限)'' を用いて描画できる。
''abline(a, b)'' は $y=ax+b$ の直線を描画する関数である。
<code rsplus>
> f <- function(x) {return(x ^ 3 - 9 * x)}  # 関数f(x)を定義
> plot(f, -5, 5) # y=f(x)のグラフを描画
> abline(0, 0) # y=0 の線を描画
</code>
$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'' のループを途中で抜けるためのコマンドである。

<file rsplus bisection.R[enable_line_numbers="true"]>
# 二分法

# 関数の定義
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
  }
}
</file>

2点の初期値を $x_{1}=1$, $x_{2}=10$ としているので、これを実行すると解の一つ $x=3$ に収束すると思われるが、実際はどうだろうか。

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

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

関数（5行目）の書き換え
<code rsplus>
  return(x ^ 2 - 2)
</code>

$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}$ の値を求めることと同じである。


===== 小テスト① =====
[[KMS>|Moodle Server（非公式）]]で第3回の小テスト①を受験しなさい。

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

代数方程式 $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)}$ が発散してしまい収束が悪くなる。

{{ newton-method.png?nolink |ニュートン法}}

===== 小テスト② =====

[[KMS>|Moodle Server（非公式）]]で第3回の小テスト②を受験しなさい。



===== ニュートン法の例 =====

==== 例1 ====

関数 
\[
f(x)=x^{3}-9x
\]
を考える。

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

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

<file rsplus newton.R[enable_line_numbers="true"]>
# ニュートン法

# 関数の定義
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
  }
}
</file>


$x^{(0)}$ の初期値を $x^{(0)}=3.6$  にすると解の一つ $x=3$ に収束する。
$x^{(0)}$ の初期値を変えることで他の解に収束する。

----

==== 例2 ====

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

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

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

関数（5行目）の書き換え
<code rsplus>
  return(x ^ 3 - 9 * x + 28)
</code>

$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行目）の書き換え
<code rsplus>
  return(x ^ 2 - 4 * x + 4)
</code>
微分（10行目）の書き換え
<code rsplus>
  return(2 * x - 4)
</code>


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



===== 課題 =====
[[KMS>|Moodle Server（非公式）]] で第3回の課題を行いなさい。

締め切り：2025年4月30日（水）20時