前回紹介したニュートン・ラフソン法を利用して、方程式の解を求めてみましょう。練習問題ですから、シンプルで、手でも計算が可能な方程式を取り上げます。問題の方程式を因数分解をするとわかりますが、解は重解で1つのみです。数値計算して得た結果と、因数分解して得た結果がどれだけ一致しているか、それを確認してみましょう。
問題 ニュートン・ラフソン法で4x2+12x+9=0 の解を求めましょう。
前回の例題で紹介したソースコードを編集して、問題の方程式のためのプログラムを作成しましょう。変更が必要なのは、作成した漸化式を表すコードの部分だけです。漸化式さえコードにしてしまえば終わりの簡単な課題ではつまりませんので、もうひと味。方程式の真の値を手で計算しておいて、真の値とニュートン・ラフソン法で得た値との差を最後に表示してみましょう。
解説
問題 ニュートン・ラフソン法で4x2+12x+9=0 の解を求めましょう。
先ず、漸化式を導きます。
方程式を手で解いた結果は、x=-1.5 です。
漸化式をdouble型の精度いっぱいを目指して計算させるプログラムを次に示します。値が振動する、あるいは収束せずに発散する場合もありますので、繰り返し回数は最大で1000回に設定しました。
実行結果は次のようになりました。
解析的に方程式の解を求めると、重解でx=-1.5 です。今回問題で取り上げた方程式の解の精度は小数点以下8桁までといったところです。残念ながら、これ以上の精度は得られませんでした。double型の精度ぎりぎりまで正確な値が得られることを期待したいところですが、浮動小数点数の誤差の影響でこの程度が限界となっていました。
少々悔しいので精度を改善できないか検討してみましょう。式68.2をよく見ると、xn+1はxn から変化分を「削って」作っていることがわかります。この削る部分が小さくなりすぎると、プログラムは計算を終了します。
考えられる誤差の発生機構は、減算の時に発生する情報欠落(連載第12回)のようです。ごく近い値の減算による桁落ち(連載第10回)ではありません。
試しに、桁落ちによる誤差を防ぐ操作を施してみます。
aはdouble型の一時変数です。このように計算式を変形して実行してみた結果は次の通りでした。
誤差はわずかに大きくなってしまいました。計算回数が増えたことにより、丸め誤差(連載第12回)が多く積み重なってしまったのでしょう。正しい改善策を選ばないと、かえって悪い結果を得るという良い例です。
では、丸め誤差を防ぐにはどうすればよいでしょうか。これ以上精度の良い変数は使えない、つまり変数の精度はこれ以上上げられないとすると、丸め操作の回数を減らすことが最善の策のようです。式変形をして、計算操作をなるべく少なくしてみましょう。
漸化式をこのように式変形することが出来ました。こうしたときの実行結果は次の通りでした。
やりました! double型の精度いっぱいまで追い込むことに成功しました。真の値が分かっている場合には、このように自信を持って式変形と検証に取り組むことが出来ます。しかし、そもそも真の値が分からない場合にニュートン・ラフソン法を用いるのですから、いつもこううまくいくものではありません。
結局は、パラメーターや計算式を操作し、計算過程を画面に逐次表示し、目で確認しながら精度を上げていくのがニュートン・ラフソン法のベターな用い方でしょう。
今回はここまで
ニュートン・ラフソン法を利用して方程式を解いてみました。今回の問題に取り組んでみて「意外と簡単だった」と感じていただければ幸いです。より高次の方程式の解を求めるときには、初期値の決め方や終了条件、誤差への対応など、考えることが多くあります。そもそも、ニュートン・ラフソン法が向かない方程式だって多くあります。本格的な活用が必要になったら専門書にあたり、先人達が積み上げたノウハウを有り難く利用しましょう。自分でゼロから取り組むのはとても大変ですし、時間がもったいないことです。