gihyo.jpサイトのロゴ

はじめMath! Javaでコンピュータ数学

第68回微分・積分の数学　ニュートン・ラフソン法　［後編］

2009-10-08

前回紹介したニュートン・ラフソン法を利用して、方程式の解を求めてみましょう。練習問題ですから、シンプルで、手でも計算が可能な方程式を取り上げます。問題の方程式を因数分解をするとわかりますが、解は重解で1つのみです。数値計算して得た結果と、因数分解して得た結果がどれだけ一致しているか、それを確認してみましょう。

問題　ニュートン・ラフソン法で4x2＋12x＋9＝0 の解を求めましょう。

前回の例題で紹介したソースコードを編集して、問題の方程式のためのプログラムを作成しましょう。変更が必要なのは、作成した漸化式を表すコードの部分だけです。漸化式さえコードにしてしまえば終わりの簡単な課題ではつまりませんので、もうひと味。方程式の真の値を手で計算しておいて、真の値とニュートン・ラフソン法で得た値との差を最後に表示してみましょう。

解説

問題　ニュートン・ラフソン法で4x2＋12x＋9＝0 の解を求めましょう。

先ず、漸化式を導きます。

方程式を手で解いた結果は、x＝－1.5 です。

漸化式をdouble型の精度いっぱいを目指して計算させるプログラムを次に示します。値が振動する、あるいは収束せずに発散する場合もありますので、繰り返し回数は最大で1000回に設定しました。

ソースコード：NewtonMethodTest2.java
//filename : NewtonMethodTest2.java

public class NewtonMethodTest2 {

  //---------------------------------
  //データ宣言部
  static int ShokiTi=1;
  static int Kurikaeshi=1000;

  /**
   * メインメソッド
   */
  public static void main (String args []) {
    NewtonMethod();
  }// end of main()


  //---------------------------------
  //メイン関数内で利用する関数宣言
  //- - - - - - - - - - - - - - - - -
  //
  static void NewtonMethod(){
    double x_p;//漸化式の右側のx 値
    double x_s;//漸化式の左側のx値
    x_p = ShokiTi;
    for(int i=0;i<=Kurikaeshi;i++){
      x_s = x_p - (4*x_p*x_p+12*x_p+9) / (8*x_p+12);
      System.out.printf
         ("%4d: %20.18f -> %20.18f\n",i, x_p , x_s );
      if (x_s == x_p) {
        System.out.println
           ("精度の限界まで計算しました。e="+ Math.abs(x_s+1.5));
        break;
      }// of if
      x_p = x_s;
    }// of for i

  }// end of NewtonMethod

}// end of this file...NewtonMethodTest2.java

実行結果は次のようになりました。

画面
C:\>java NewtonMethodTest2
   0: 1.000000000000000000 -> -0.250000000000000000
   1: -0.250000000000000000 -> -0.875000000000000000
   2: -0.875000000000000000 -> -1.187500000000000000
   3: -1.187500000000000000 -> -1.343750000000000000
   4: -1.343750000000000000 -> -1.421875000000000000
   5: -1.421875000000000000 -> -1.460937500000000000
   6: -1.460937500000000000 -> -1.480468750000000000
   7: -1.480468750000000000 -> -1.490234375000000000
   8: -1.490234375000000000 -> -1.495117187500000000
   9: -1.495117187500000000 -> -1.497558593750000000
   10: -1.497558593750000000 -> -1.498779296875000000
   11: -1.498779296875000000 -> -1.499389648437500000
   12: -1.499389648437500000 -> -1.499694824218750000
   13: -1.499694824218750000 -> -1.499847412109375000
   14: -1.499847412109375000 -> -1.499923706054687500
   15: -1.499923706054687500 -> -1.499961853027343800
   16: -1.499961853027343800 -> -1.499980926513671900
   17: -1.499980926513671900 -> -1.499990463256836000
   18: -1.499990463256836000 -> -1.499995231628418000
   19: -1.499995231628418000 -> -1.499997615814209000
   20: -1.499997615814209000 -> -1.499998807907104500
   21: -1.499998807907104500 -> -1.499999403953552200
   22: -1.499999403953552200 -> -1.499999701976776100
   23: -1.499999701976776100 -> -1.499999850988388000
   24: -1.499999850988388000 -> -1.499999925494194000
   25: -1.499999925494194000 -> -1.499999961256981000
   26: -1.499999961256981000 -> -1.499999972719412700
   27: -1.499999972719412700 -> -1.499999980858702300
   28: -1.499999980858702300 -> -1.499999992458992400
   29: -1.499999992458992400 -> -1.499999992458992400
精度の限界まで計算しました。e=7.541007596145732E-9

解析的に方程式の解を求めると、重解でx＝－1.5 です。今回問題で取り上げた方程式の解の精度は小数点以下8桁までといったところです。残念ながら、これ以上の精度は得られませんでした。double型の精度ぎりぎりまで正確な値が得られることを期待したいところですが、浮動小数点数の誤差の影響でこの程度が限界となっていました。 少々悔しいので精度を改善できないか検討してみましょう。式68.2をよく見ると、xn+1はxn から変化分を「削って」作っていることがわかります。この削る部分が小さくなりすぎると、プログラムは計算を終了します。

考えられる誤差の発生機構は、減算の時に発生する情報欠落（連載第12回）のようです。ごく近い値の減算による桁落ち（連載第10回）ではありません。

試しに、桁落ちによる誤差を防ぐ操作を施してみます。

ソースコード：NewtonMethodTest3.java部分
a = (4*x_p*x_p+12*x_p+9) / (8*x_p+12);
x_s = (x_p*x_p - a*a)/(x_p+a);

aはdouble型の一時変数です。このように計算式を変形して実行してみた結果は次の通りでした。

画面
    (略)
  40: -1.500000008809426000 -> -1.500000008809426000
精度の限界まで計算しました。e=8.809426077505123E-9

誤差はわずかに大きくなってしまいました。計算回数が増えたことにより、丸め誤差（連載第12回）が多く積み重なってしまったのでしょう。正しい改善策を選ばないと、かえって悪い結果を得るという良い例です。

では、丸め誤差を防ぐにはどうすればよいでしょうか。これ以上精度の良い変数は使えない、つまり変数の精度はこれ以上上げられないとすると、丸め操作の回数を減らすことが最善の策のようです。式変形をして、計算操作をなるべく少なくしてみましょう。

ソースコード：NewtonMethodTest4.java部分
x_s = x_p - (0.5*x_p+3.0/4.0);

漸化式をこのように式変形することが出来ました。こうしたときの実行結果は次の通りでした。

画面
    (略)
  54: -1.500000000000000000 -> -1.500000000000000000
精度の限界まで計算しました。e=0.0

やりました！ double型の精度いっぱいまで追い込むことに成功しました。真の値が分かっている場合には、このように自信を持って式変形と検証に取り組むことが出来ます。しかし、そもそも真の値が分からない場合にニュートン・ラフソン法を用いるのですから、いつもこううまくいくものではありません。

結局は、パラメーターや計算式を操作し、計算過程を画面に逐次表示し、目で確認しながら精度を上げていくのがニュートン・ラフソン法のベターな用い方でしょう。

今回はここまで

ニュートン・ラフソン法を利用して方程式を解いてみました。今回の問題に取り組んでみて「意外と簡単だった」と感じていただければ幸いです。より高次の方程式の解を求めるときには、初期値の決め方や終了条件、誤差への対応など、考えることが多くあります。そもそも、ニュートン・ラフソン法が向かない方程式だって多くあります。本格的な活用が必要になったら専門書にあたり、先人達が積み上げたノウハウを有り難く利用しましょう。自分でゼロから取り組むのはとても大変ですし、時間がもったいないことです。

おすすめ記事

記事・ニュース一覧