ゆるゆるプログラミング

・ニュートン法で平方根 その2

ニュートン法は、関数f(x)が与えられたときf(x)=0となるxを数値計算の反復によって求めるアルゴリズムの1つです。その計算式は以下です。

ニュートン法 漸化式

ここでは、数値微分を使ったニュートン法平方根を計算するプログラムを紹介します。

ニュートン法の式の中にf(x)をxで微分したf'(x)を含んでいますが、関数f(x)によってf'(x)を導くことが困難または不可能なことがあります。

そこで、微分の関数が分からなくてもf'(x)が計算できる数値微分を使い、ニュートン法の計算を行えるようにしました。

関数をf(x)=x2-aとし、f(x)=0となるxを求めることで、aの平方根を計算します。

この式は、以下の手順で導いています。

ニュートン法 平方根を求める元の式

以下がそのJavaソースコードです。

NewtonsMethod2.java ← クリックしてダウンロードページに移動
001:    public class NewtonsMethod2 {
002:    	// f(x)=x^2-a
003:    	static double f( double x, double a )
004:    	{
005:    		return x * x - a;
006:    	}
007:    
008:    
009:    	// f(x)の微分f'(x) - 接線の傾き
010:    	// 微分の式が分からない場合、微分の近似値を計算
011:    	static double df_ex( double x, double a )
012:    	{
013:    		double h = 0.0000001;
014:    
015:    		return ( f( x + h, a ) - f( x, a ) ) / h;
016:    	}
017:    
018:    
019:    	// ニュートン法 漸化式
020:    	static double newtons( double x, double a )
021:    	{
022:    		return x - f( x, a ) / df_ex( x, a );
023:    	}
024:    
025:    
026:    	// メイン
027:    	public static void main(String[] args) {
028:    		// 変数の宣言
029:    		double e = 0.000001;	// 収束の判定値
030:    		double x;		// 計算値を格納
031:    		double x_ans;		// 1つ前の計算値
032:    		double a;		// 平行根を求める値
033:    
034:    		// 平行根を求める値
035:    		a = 169.0;	// 13×13
036:    
037:    		// 初期値を代入
038:    		x_ans = 10000.0;
039:    
040:    		// 無限ループ
041:    		for ( int i = 1; ; ++ i ) {
042:    			// ニュートン法で計算
043:    			x = newtons( x_ans, a );
044:    
045:    			// 収束の判定
046:    			if ( e > Math.abs( x - x_ans ) ) break;
047:    
048:    			// 最新の計算結果を次の計算結果に代入
049:    			x_ans = x;
050:    
051:    			// 途中経過の表示
052:    			System.out.println( "回数:" + i + "  計算結果: " + x_ans );
053:    		}
054:    
055:    		// 改行
056:    		System.out.println();
057:    
058:    		// 結果の表示
059:    		System.out.println( "計算結果: " + x_ans );
060:    
061:    		// 結果の絶対値
062:    		x_ans = Math.abs( x_ans );
063:    		System.out.println( a + "の絶対値: " + x_ans );
064:    
065:    		// f(x)=0になるかを確認
066:    		System.out.println( "f(" + x_ans + ")=" + f( x_ans, a ) );
067:    	}
068:    }

NewtonsMethod2を実行

C:\talavax\javasample>java NewtonsMethod2

169の平方根ニュートン法で計算してます。計算の途中経過と結果を表示しています。

NewtonsMethod2.javaの出力結果

回数:1  計算結果: 5000.055835163439
回数:2  計算結果: 2500.0405926789967
回数:3  計算結果: 1250.05263083155
回数:4  計算結果: 625.0944750620009
回数:5  計算結果: 312.682233554483
回数:6  計算結果: 156.6113993198339
回数:7  計算結果: 78.84524510893075
回数:8  計算結果: 40.494340546011806
回数:9  計算結果: 22.333882320006264
回数:10  計算結果: 14.950430379777575
回数:11  計算結果: 13.127226391861301
回数:12  計算結果: 13.000616525291765
回数:13  計算結果: 13.000000014619735

計算結果: 13.000000014619735
169.0の絶対値: 13.000000014619735
f(13.000000014619735)=3.8011310721230984E-7

001:    public class NewtonsMethod2 {

クラス名を、NewtonsMethod2としています。

002:    	// f(x)=x^2-a
003:    	static double f( double x, double a )
004:    	{
005:    		return x * x - a;
006:    	}

ニュートン法でf(x)=0を求めたい関数です。ここでは、f(x)=x2-aという関数を記述しています。平方根の元の値をaに代入します。

009:    	// f(x)の微分f'(x) - 接線の傾き
010:    	// 微分の式が分からない場合、微分の近似値を計算
011:    	static double df_ex( double x, double a )
012:    	{
013:    		double h = 0.0000001;
014:    
015:    		return ( f( x + h, a ) - f( x, a ) ) / h;
016:    	}

f(x)の微分式f'(x)を数値微分で行うメソッドです。xが0.0000001増加したときの、yの増加分を計算しています。f'(x)は(yの増加分/xの増加分)なので、その値をreturn文で戻しています。ここでは、xの増加分を0.0000001としていますが関数f(x)によって変更する必要があります。いろいろな値で試してください。

019:    	// ニュートン法 漸化式
020:    	static double newtons( double x, double a )
021:    	{
022:    		return x - f( x, a ) / df_ex( x, a );
023:    	}

ニュートン法の漸化式のメソッドです。

026:    	// メイン
027:    	public static void main(String[] args) {

このmainメソッドからプログラムを実行します。

028:    		// 変数の宣言
029:    		double e = 0.000001;	// 収束の判定値
030:    		double x;		// 計算値を格納
031:    		double x_ans;		// 1つ前の計算値
032:    		double a;		// 平行根を求める値

このプログラムで使う変数を宣言しています。double型変数eは、漸化式を終了させるための値で、計算したxn+1とxn絶対値の差がe未満になったときに計算を終了させます。

034:    		// 平行根を求める値
035:    		a = 169.0;	// 13×13

平方根を求める値をdouble型変数aに代入しています。

ここでは169.0を変数aに代入しています。これにより、ニュートン法で収束させた値が169.0の平方根である13.0になることが期待できます。

037:    		// 初期値を代入
038:    		x_ans = 10000.0;

計算の初期値に10000.0を代入しています。この値は任意です。

040:    		// 無限ループ
041:    		for ( int i = 1; ; ++ i ) {

for文無限ループを作っています。int型変数iを1からインクリメント(1を足す)しています、終了条件を記述していないので無限ループになります。

042:    			// ニュートン法で計算
043:    			x = newtons( x_ans, a );

ニュートン法の漸化式にx_ansとaを渡して、次のxの値を取得しています。

045:    			// 収束の判定
046:    			if ( e > Math.abs( x - x_ans ) ) break;

計算した値xと、計算に使った値x_ansの差の絶対値を計算し、その値がeより小さければ計算が収束したと判断し、break文for文を抜けます。

048:    			// 最新の計算結果を次の計算結果に代入
049:    			x_ans = x;

x_ansにxを代入しています。

055:    		// 改行
056:    		System.out.println();
057:    
058:    		// 結果の表示
059:    		System.out.println( "計算結果: " + x_ans );

ニュートン法で求めた結果x_ansをコンソールに出力しています。

061:    		// 結果の絶対値
062:    		x_ans = Math.abs( x_ans );
063:    		System.out.println( a + "の絶対値: " + x_ans );

計算結果x_ansの絶対値をとった値をコンソールに出力しています。これが、実際に求めたい平方根の値になります。このプログラムの場合、f(x)=x2-169.0を0.0にするxを求めるようにしてるため、xの値が13.0と-13.0の2つの値のどちらでもf(x)=0を満たすことが出来ます。これは初期値によりプラスとマイナスのどちらかの値をとります。例えば初期値を-10000.0にするとx_ans=-13になります。

065:    		// f(x)=0になるかを確認
066:    		System.out.println( "f(" + x_ans + ")=" + f( x_ans, a ) );

計算結果x_ansをメソッドfに代入して、その結果が0になっているか確認するための出力です。この例では、結果の表示が指数表記になっていて、数値の後に'E-7'が付いています。これは10の-7乗を意味しています。この結果から求めたx_ansでメソッドfが戻す値は0に近く、x_ansはaの平方根であることが確認できます。

■関連コンテンツ

ニュートン法 ニュートン法について
ニュートン法で平方根 ニュートン法で平方根を計算

■新着情報

2019.09.13 長さの単位変換 1マイル、1フィートは何m?
2019.09.06 クイックソート 高速に配列に並び替える方法
2019.09.05 中央値(メディアン) 配列に格納されている値の中央値を求める
2019.09.05 最頻値 配列から出現回数が最も多い値の取得
2019.09.03 配列値の反転 配列の反転処理
2019.08.05 トランプの操作 トランプを操作するクラス

■広告

法人向けのETC専用カード

~約8,000名の受講生と80社以上の導入実績~ 企業向けプログラミング研修ならCodeCamp

日本最大級ショッピングサイト!お買い物なら楽天市場

Topへ