前回の授業では二分法のアルゴリズムについてまとめた上でプログラムとしての実装まで行いましたが、実はそのアルゴリズムには一つの欠陥があります。
入力:関数$f$、区間の両端 $a$ と $b$、許容誤差 $TOL$
出力:反復回数 $n$、近似解 $p$
Step 1:関数 $f$ を定義する。$a$、$b$、$TOL$ を初期化する。$n=0$ とする。
Step 2:$b-a\geq 2\cdot TOL$ の間、次のStep 2-1~2-3を繰り返す。
Step 2-1:$n$ に $1$ を加算する。
Step 2-2:$p=\frac{a+b}{2}$ とする。
Step 2-3:もし $f(a)\cdot f(p)<0$ なら、$b=p$ とする。そうでなければ、$a=p$ とする。
Step 3:$p=\frac{a+b}{2}$ とする。
Step 4:反復回数 $n$ と近似解 $p$ を出力する。
上記が前回説明したアルゴリズムです。
Step 2-3でもし $f(p)==0$ の場合には $a=p$ として区間の左端を更新することになり、本来は次の反復で区間の右端を更新する必要がありますが、$f(a)==0$ のため次の反復でも区間の左端を更新することになります。結果として、区間の中に厳密解の存在が保証された状態を維持できなくなります。
ここで、プログラム上で $f(p)==0$ の条件判定を行うことを考えると、小数(浮動小数点数)の計算では多くの場合に誤差が含まれるので、そのまま等価演算子==
を用いて判定することはあまり意味がありません。代わりに「ほぼ等しい」かどうかを判定するために、非常に小さな正の数 $\varepsilon$ に対して
を使うことが一般的です。
以上を踏まえると、二分法のアルゴリズムは次のように修正されます(赤字は内容の追加、青字は順序の変更)。
入力:関数$f$、区間の両端 $a$ と $b$、許容誤差 $TOL$、微小正数 $EPS$(epsilonの略)
出力:反復回数 $n$、近似解 $p$
Step 1:関数 $f$ を定義する。$a$、$b$、$TOL$、$EPS$ を初期化する。$n=0$ とする。
Step 2:$p=\frac{a+b}{2}$ とする。
Step 3:$b-a\geq 2\cdot TOL$ の間、次のStep 3-1~3-3を繰り返す。
Step 3-1:もし $|f(p)|
Step 3-2:$n$ に $1$ を加算する。
Step 3-3:$p=\frac{a+b}{2}$ とする。
Step 4:反復回数 $n$ と近似解 $p$ を出力する。
改めて前回の演習2のコードを書くと、次のようになります。ただし、微小正数 $EPS=10^{-15}$ とします。
function value = f(x)
value = x.^3-2;
end
a = 0;
b = 3;
TOL = 0.00001; %下の行と同様に1E-5と書いてもよい
EPS = 1E-15; %1E-15は指数表記であり、1*10^(-15)と同じ意味
n = 0;
p = (a+b)/2;
while b-a >= 2*TOL
if abs(f(p)) < EPS %f(p)が0とほぼ等しいかどうかを判定する absは絶対値を求める関数
break %while文による繰り返しを強制終了する
elseif f(a)*f(p) < 0
b = p;
else
a = p;
end
n = n+1;
p = (a+b)/2;
end
printf("反復回数n=%d, 近似解p=%f\n",n,p)
反復回数n=18, 近似解p=1.259920
指数表記(exponential notation)は、科学技術分野でよく使用される $a\times 10^b$ の形の表記法であり、非常に大きな数や $0$ に近い数を表現する際に便利です。
MATLABを含むコンピュータ上では、仮数部 $a$ と指数部 $b$ の間にE
またはe
を挟んで表記します。
1E-15
1e-15 %大文字でも小文字でも同じ
1.234E+5
1.234E5 %指数部の+は省略できる
1.234E-5
ans = 1.0000e-15 ans = 1.0000e-15 ans = 123400 ans = 123400 ans = 0.000012340
ニュートン法(Newton's method)は、二分法と同様に方程式の近似解を反復計算によって求める有名なアルゴリズムです。ニュートン・ラフソン法(Newton-Raphson method)と呼ばれることもあります。
微分可能な関数 $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$ の近似解として得ることができます。
ニュートン法における反復の終了条件はいくつか考えられますが、生成した数列の項がほぼ変化しなくなったとき、すなわち、微小正数 $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}$ の役割として使用すればそれで十分です。
以上から、ニュートン法のアルゴリズムは次のようにまとめられます。
入力:関数$f$、$f$ の微分 $f\_dif$、初期値 $x\_old=x_0$、微小正数 $EPS$、反復回数の上限 $N$
出力:反復回数 $n$ と近似解 $x\_new$、または「反復回数の上限を超えた」旨のメッセージ
Step 1:関数 $f$ と $f\_dif$ を定義する。$x\_old$、$EPS$、$N$ を初期化する。
Step 2:$x\_old$、$f$、$f\_dif$ を用いて $x\_new$ を初期化する。$n=1$ とする。
Step 3:$|x\_new-x\_old|\geq EPS$ かつ $n\leq N$ の間、次のStep 3-1~3-3を繰り返す。
Step 3-1:$x\_old=x\_new$ とする。
Step 3-2:$x\_old$、$f$、$f\_dif$ を用いて $x\_new$ を更新する。
Step 3-3:$n$ に $1$ を加算する。
Step 4:もし $n>N$ なら、「反復回数の上限を超えた」旨のメッセージを出力する。そうでなければ、反復回数 $n$ と近似解 $x\_new$ を出力する。
赤字の箇所は敢えて少し曖昧に書いていますので、どうしたらよいか考えてください。
上記のニュートン法のアルゴリズムを実現するコードを書き、関数 $f(x)=x^2-2$ に対して方程式 $f(x)=0$ の近似解を求めてください。ただし、初期値 $x_0=0.5$、微小正数 $EPS=10^{-10}$、反復回数の上限 $N=20$ とします。
%演習1のコード
function value = f(x)
value = x^2-2;
end
function value = f_dif(x)
value = 2*x;
end
%ここに、不足しているコードを書いてください
if n > N
printf("%d回の反復では、近似解を求めることができなかった。\n",N)
else
printf("%d回の反復で、近似解%.15fを求めることができた。\n",n,x_new)
end
二分法とニュートン法はともに方程式の近似解を反復計算によって求める有名なアルゴリズムですが、比較すると次のような特徴があります。
二分法は、ニュートン法と比べて収束速度が遅いです(一次収束)。一方で、連続関数という緩い条件の下で関数値が異符号となるように区間の両端を取りさえすれば、確実に厳密解に収束するという長所があります。
ニュートン法は、収束速度が非常に速いです(二次収束)。収束速度に関する詳しい議論は「数値解析A」の授業でされているため、ここでは控えます。一方で、与えられた関数の微分を求める必要があり、さらに与えられた関数と初期値によって収束しない場合があることは短所と言えます。なお、ニュートン法における数列が収束するための十分条件を与える数学的な定理がいくつか知られていますが、ここでは扱いません。
ニュートン法が収束しないような関数と初期値を与えて、それを確認するコードを書いてください。特に、生成された数列の各項($x\_new$)を画面に出力して、収束しない様子が分かるようにしてください。
%演習2のコード
前述のように、ニュートン法では与えられた関数の微分を求める必要があります。関数が微分可能でない場合や関数の微分を求めることが難しい場合にも適用できるようにニュートン法を修正したアルゴリズムとして、セカント法(secant method)があります。割線法とも呼ばれています。
関数 $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$ の近似解として得ることができます。
セカント法は、ニュートン法における微分 $f'(x_{n-1})$ の代わりに、近似的に差分商(変化の割合)
$$ \frac{f(x_{n-1})-f(x_{n-2})}{x_{n-1}-x_{n-2}} $$を用いていると解釈できます。
幾何学的には、曲線の接線ではなく、代替的に割線(曲線と二点以上で交わる直線のことを割線と呼ぶ)を用いた方法と言えます。
セカント法のアルゴリズムは、ニュートン法と共通している部分が多いです。反復の終了条件は全く同じで、微小正数 $EPS$ と反復回数の上限 $N$ に対して
$$ |x_n-x_{n-1}|<EPS $$または
$$ n>N $$が成り立つこととします。
ニュートン法と違う点として、三つの変数 $x\_new$、$x\_old$、$x\_old2$ を用意してそれぞれを $x_n$、$x_{n-1}$、$x_{n-2}$ の役割として使用すべきことに注意してください。
セカント法のアルゴリズムを実現するコードを書き、関数 $f(x)=x^2-2$ に対して方程式 $f(x)=0$ の近似解を求めてください。ただし、初期値 $x_{-1}=0.3$、$x_0=0.5$、微小正数 $EPS=10^{-10}$、反復回数の上限 $N=20$ とします。
%演習3のコード
演習1~演習3に取り組んでください。