プログラミング演習第10回方程式を解く（二分法，ニュートン法）

目次

1 二分法 1 2 ニュートン法 2 3 実習課題 3 4 練習問題 4

図 1:

1 二分法

方程式 f(x) = 0 の実数解は，関数 y = f(x) のグラフと， x 軸との共有点の x 座標にほかならない．

関数のグラフを調べて，方程式の実数解の近似値，すなわち近似解を求めることを考えよう．

一般に，関数 y ＝ f(x) のグラフが，ある区間でとぎれることなく，つながっているとき，関数 f(x) は，その区間で連続であるという．区間 [a, b] で連続な関数 f ＝ (x) について，次の定理が成り立つ．

定理 f(a) と f(b) の値の符号が異なるとき，すなわち， f(a)f(b) < 0 ならば， f(a) = 0 を満たす実数 α が， a と b の間に少なくとも１つ存在する．

この定理を利用して，方程式の近似解を求めることを考えよう．関数 f(x) は，区間 [a, b] において連続で， f(a)f(b) < 0 を満たすとし，また，方程式 f(x) ＝ 0 はこの区間でただ１つの解をもつと仮定する．

このとき，区間 [a, b] の中点を c = \frac{a+b}{2} とすると f(c) の値は 0 となるか，または f(a),f(b) のいずれかと異符号である．

ここで， f(c) = 0 なら， c は方程式の解である．

また， f(a)f(c) < 0 ならば, 方程式は区間 [a, c] に解をもち， f(b)f(c) < 0 ならば，方程式は区間 [c, b] に解をもつ．

このように，区間の中点を考えることにより，解の存在する範囲を初めの半分に縮小することができる．そして，この操作を繰り返すと，誤差の限界を必要なだけ小さくした近似解が得られる．

このようにして近似解を求める方法を二分法という．

一般に，方程式 f(x) = 0 の近似解を，誤差の限界は ϵ として二分法で求める手順は次のようになる．

(1) 変数 a, b に数値を入力して， f(a)f(b) < 0 を満たすかどうか判断し，正しければ (2) にすすみ，正しくなければ計算不能と表示する．

図 2:

(2) 区間 [a, b] の中点 c = \frac{a+b}{2} を求める．

(3) f(c) = 0 ならば， c を解とする． f(a)f(c) < 0 ならば，変数 b に c の値を代入し， f(a)f(c) > 0 ならば，変数 a に c の値を代入する．

(4) |a − b| ≤ ε かどうかを判断し，正しければ c の値を近似解とする．正しくなければ (2) に戻り，操作を続ける．

上の (3) において， f(c) ≠ 0 のとき，変数 a, b の値をこのようにおき換えると， f(a)f(b) < 0 となり，真の解 α は区間 [a, b] にある．また，このとき， c の値は a か b に等しいから， |c − α| ≤ |a − b| となり， c を近似解としたときの誤差の限界は |a − b| である．したがって， (4) で誤差の限界が判定される．

方程式 f(x) = 0 が区間 [a, b] で 3 つの実数解 l, m, n (a < l < m < c < n < b) をもち, f(a) > 0, f(c) > 0, f(b) < 0 とする． c を a と b の中点として二分法を適用すると，第２段階以後区間 [a, c] は無視されて，区間 [c, b] にある解だけを近似することになる．このように，ある区間で方程式が２つ以上の解をもつ場合，二分法で近似解を求めると， 1 つの解以外は無視される結果になる．また，例えば，方程式 (x − 1)^{2} = 0 の重解１の近似値は，二分法では求められない．したがって，二分法を適用する場合，あらかじめ関数のグラフの概形を書くなどして，解について調べておくことが必要である．

2 ニュートン法

方程式の近似解を求める方法には，前項で学んだ二分法のほかに，曲線の接線を利用する方法がある．

任意の実数 a に対して，関数 f(x) の x = a における微分係数 f′(a) が存在する場合，方程式 f(x) = 0 の解 α の近似値を求めることを考えよう．

そのため，まず α に近い値 x_{1} を１つとり，関数 y = f(x) のグラフ上の点 (x_{1}, f(x_{1})) における接線が， x 軸と交わる点の x 座標を x_{2} とする．

このとき，接線は，点 (x_{1}, f(x_{1})) を通り，傾きが f′(x_{1}) の直線であるから，その方程式は次のようになる．

y = f′(x_{1})(x − x_{1}) + f(x_{1})

ここで， y = 0 とおいて， x について解くと x_{2} が求まる．すなわち

x_{2} = x_{1} − \frac{f(x_{1})}{f′(x_{1})}  (1)

次に，点 (x_{2}, f(x_{2})) における曲線の接線が， x 軸と交わる点の x 座標を x_{3} とし，以下この操作を繰り返すと，数列

x_{1}, x_{2}, x_{3}, ⋯, x_{n}, ⋯  (2)

が得られる．そして， (1) と同様にして， n = 2, 3, ⋯ に対しても等式

x_{n+1} = x_{n} − \frac{f(x_{n})}{f′(x_{n})}  (3)

が成り立ち， (2) は漸化式 (3) で定まる．

この数列 (2) は，最初の値 x_{1} が α に十分近いとき，多くの場合急速に α に近づき，これを利用して方程式の近似解が効率よく求められる．

このようにして近似解を求める方法をニュートン法といい，最初にとった値 x_{1} を初期値という．

ニュートン法では，二分法と違って，誤差の限界の正確な判定が難しい．そのため，数列 (2) において， |x_{n} − x_{n−1}| < ϵ となったとき， x_{n} は誤差の限界が ϵ の近似解であると考える．

方程式 f(x) = 0 の近似解を，誤差の限界は ϵ として，ニュートン法で求める手順をまとめると，次のようになる．

(1) 適当な初期値 a を入力する．

(2) b = a − \frac{f(a)}{f′(a)} を求める．

(3) |b − a| < ϵ かどうかを判断し，正しければ， b を近似解とする．正しくなければ，変数 a に b の値を入力して， (2) にもどる．

3 実習課題

以下の問題を解き, 解答を時間内に大坂にメールで送ること. 但し, Subject名は, creportno.10-1 とする.

時間内にできない場合は 6 月 21 日までに提出せよ.

(1) 初期値を 0, 2(即ち, 区間 [0,2]), 誤差を限界 ϵ を ϵ = 10^{−6} = 0.000001 として, 二分法で x^{2} = 2 の解を出力する下記のプログラムを完成せよ. また, そのとき得られた解を書け. (Hint:1 ページ目の二分法の手順を参照せよ)

#include <stdio.h>
double f(double);
main()
{
  double a, b, c, e;
  a = (1);
  b = (2);
  e= 0.000001;                      /* 誤差 */
while((3))
{
  if((4))
     c = (5);
   else
    printf("計算不能");
  if(f(c)== 0)
   printf("求める解は %f です",c);
  if((6))
    b = (7);
      else if((8))
   a = (9);
}
    printf("近似解は %f です", c);
}
double f(double x)
{
   (10);
 }

4 練習問題

次の問題を解答したものは評価する. 但し, subject 名は creportno.10-2 とし, 締切りは 6 月 28 日とする.

• 初期値 k を入力し，ニュートン法で x^{3} = k の解を出力するプログラムを作成せよ．ただし誤差の限界 ϵ は ϵ = 10^{−6} = 0.000001 とする．