2022年度プログラミング演習A・B

第5回:非線形方程式の解の計算2

0. 前回の補足

前回の授業では二分法のアルゴリズムについてまとめた上でプログラムとしての実装まで行いましたが、実はそのアルゴリズムには一つの欠陥があります。



上記が前回説明したアルゴリズムです。

Step 2-3でもし $f(p)==0$ の場合には $a=p$ として区間の左端を更新することになり、本来は次の反復で区間の右端を更新する必要がありますが、$f(a)==0$ のため次の反復でも区間の左端を更新することになります。結果として、区間の中に厳密解の存在が保証された状態を維持できなくなります。

ここで、プログラム上で $f(p)==0$ の条件判定を行うことを考えると、小数(浮動小数点数)の計算では多くの場合に誤差が含まれるので、そのまま等価演算子==を用いて判定することはあまり意味がありません。代わりに「ほぼ等しい」かどうかを判定するために、非常に小さな正の数 $\varepsilon$ に対して

$$ |f(p)|<\varepsilon $$

を使うことが一般的です。

以上を踏まえると、二分法のアルゴリズムは次のように修正されます(赤字は内容の追加、青字は順序の変更)。



改めて前回の演習2のコードを書くと、次のようになります。ただし、微小正数 $EPS=10^{-15}$ とします。

(補足)指数表記について

指数表記(exponential notation)は、科学技術分野でよく使用される $a\times 10^b$ の形の表記法であり、非常に大きな数や $0$ に近い数を表現する際に便利です。

MATLABを含むコンピュータ上では、仮数部 $a$ と指数部 $b$ の間にEまたはeを挟んで表記します。

1. ニュートン法

ニュートン法(Newton's method)は、二分法と同様に方程式の近似解を反復計算によって求める有名なアルゴリズムです。ニュートン・ラフソン法(Newton-Raphson method)と呼ばれることもあります。

1.1. ニュートン法の考え方

微分可能な関数 $f:\mathbb{R}\to\mathbb{R}$ が与えられているとし、方程式 $f(x)=0$ の近似解を求めることを考えます。

初期値として、厳密解に近いと思われる値 $x_0$ を一つ取ります。

$x$-$y$ 平面上で曲線 $y=f(x)$ の点 $(x_0,f(x_0))$ における接線を考え、その $x$ 切片($x$ 軸との交点の $x$ 座標)を $x_1$ とします。ここで、考えている接線の方程式は

$$ y=f'(x_0)(x-x_0)+f(x_0) $$

であるため、点 $(x_1,0)$ を代入することで

$$ x_1=x_0-\frac{f(x_0)}{f'(x_0)} $$

が得られます。多くの場合に、$x_1$ は $x_0$ よりも厳密解に近い値となります。

この操作を繰り返して、漸化式

$$ x_n=x_{n-1}-\frac{f(x_{n-1})}{f'(x_{n-1})} $$

によって定まる数列 $\{x_n\}$ を生成します。多くの場合に数列 $\{x_n\}$ は厳密解に収束するため、反復を適当なところで終了すれば最後に求めた項を方程式 $f(x)=0$ の近似解として得ることができます。

1.2. プログラムとしての実装

ニュートン法における反復の終了条件はいくつか考えられますが、生成した数列の項がほぼ変化しなくなったとき、すなわち、微小正数 $EPS$ に対して

$$ |x_n-x_{n-1}|<EPS $$

を満たすときに、「数列が収束した」とみなして反復を終了することにします。ただし、近似解と厳密解の誤差が $EPS$ より小さくなることは保証されていませんので注意してください。

加えて、関数 $f$ と初期値 $x_0$ によっては数列 $\{x_n\}$ が収束しない場合もあるため、反復回数が事前に決めた上限 $N$ を超えたときには「数列が収束しなかった」とみなして反復を終了することにします。

ニュートン法のプログラムでは数列 $\{x_n\}$ の各項を配列に格納していくこともできますが、二つの変数 $x\_new$ と $x\_old$ を用意してそれぞれを $x_n$ と $x_{n-1}$ の役割として使用すればそれで十分です。

以上から、ニュートン法のアルゴリズムは次のようにまとめられます。



赤字の箇所は敢えて少し曖昧に書いていますので、どうしたらよいか考えてください。

演習1

上記のニュートン法のアルゴリズムを実現するコードを書き、関数 $f(x)=x^2-2$ に対して方程式 $f(x)=0$ の近似解を求めてください。ただし、初期値 $x_0=0.5$、微小正数 $EPS=10^{-10}$、反復回数の上限 $N=20$ とします。

1.3. 二分法とニュートン法の比較

二分法とニュートン法はともに方程式の近似解を反復計算によって求める有名なアルゴリズムですが、比較すると次のような特徴があります。

二分法は、ニュートン法と比べて収束速度が遅いです(一次収束)。一方で、連続関数という緩い条件の下で関数値が異符号となるように区間の両端を取りさえすれば、確実に厳密解に収束するという長所があります。

ニュートン法は、収束速度が非常に速いです(二次収束)。収束速度に関する詳しい議論は「数値解析A」の授業でされているため、ここでは控えます。一方で、与えられた関数の微分を求める必要があり、さらに与えられた関数と初期値によって収束しない場合があることは短所と言えます。なお、ニュートン法における数列が収束するための十分条件を与える数学的な定理がいくつか知られていますが、ここでは扱いません。

演習2(オプション)

ニュートン法が収束しないような関数と初期値を与えて、それを確認するコードを書いてください。特に、生成された数列の各項($x\_new$)を画面に出力して、収束しない様子が分かるようにしてください。

2. セカント法

前述のように、ニュートン法では与えられた関数の微分を求める必要があります。関数が微分可能でない場合や関数の微分を求めることが難しい場合にも適用できるようにニュートン法を修正したアルゴリズムとして、セカント法(secant method)があります。割線法とも呼ばれています。

2.1. セカント法の考え方

関数 $f:\mathbb{R}\to\mathbb{R}$ が与えられているとし、方程式 $f(x)=0$ の近似解を求めることを考えます。

初期値として、厳密解に近いと思われる二つの異なる値 $x_{-1}$ と $x_0$ を取ります。

$x$-$y$ 平面上で点 $(x_{-1},f(x_{-1}))$ と $(x_0,f(x_0))$ を通る直線を考え、その $x$ 切片($x$ 軸との交点の $x$ 座標)を $x_1$ とします。ここで、考えている直線の方程式は

$$ y=\frac{f(x_0)-f(x_{-1})}{x_0-x_{-1}}(x-x_0)+f(x_0) $$

であるため、点 $(x_1,0)$ を代入することで

$$ x_1=x_0-\frac{f(x_0)(x_0-x_{-1})}{f(x_0)-f(x_{-1})} $$

が得られます。多くの場合に、$x_1$ は $x_{-1}$ と $x_0$ よりも厳密解に近い値となります。

この操作を繰り返して、漸化式

$$ x_n=x_{n-1}-\frac{f(x_{n-1})(x_{n-1}-x_{n-2})}{f(x_{n-1})-f(x_{n-2})} $$

によって定まる数列 $\{x_n\}$ を生成します。多くの場合に数列 $\{x_n\}$ は厳密解に収束するため、反復を適当なところで終了すれば最後に求めた項を方程式 $f(x)=0$ の近似解として得ることができます。

2.2. ニュートン法に対するセカント法の解釈

セカント法は、ニュートン法における微分 $f'(x_{n-1})$ の代わりに、近似的に差分商(変化の割合)

$$ \frac{f(x_{n-1})-f(x_{n-2})}{x_{n-1}-x_{n-2}} $$

を用いていると解釈できます。

幾何学的には、曲線の接線ではなく、代替的に割線(曲線と二点以上で交わる直線のことを割線と呼ぶ)を用いた方法と言えます。

2.3. プログラムとしての実装

セカント法のアルゴリズムは、ニュートン法と共通している部分が多いです。反復の終了条件は全く同じで、微小正数 $EPS$ と反復回数の上限 $N$ に対して

$$ |x_n-x_{n-1}|<EPS $$

または

$$ n>N $$

が成り立つこととします。

ニュートン法と違う点として、三つの変数 $x\_new$、$x\_old$、$x\_old2$ を用意してそれぞれを $x_n$、$x_{n-1}$、$x_{n-2}$ の役割として使用すべきことに注意してください。

演習3

セカント法のアルゴリズムを実現するコードを書き、関数 $f(x)=x^2-2$ に対して方程式 $f(x)=0$ の近似解を求めてください。ただし、初期値 $x_{-1}=0.3$、$x_0=0.5$、微小正数 $EPS=10^{-10}$、反復回数の上限 $N=20$ とします。

第5回レポート課題

演習1~演習3に取り組んでください。