ニュートン・ラフソン法による逆数の求め方メモ

そもそもニュートン・ラフソン法ってなんだっけ？

ニュートン・ラフソン法（単にニュートン法ともいう）とは、方程式の解を近似値として反復的に求める手法です。

ニュートン・ラフソン法で逆数や平方根などを求める場合、普通はの解がその求めたい数になるようにします。

ということで最終目標は、を満たすときのを求めることです。

大雑把な流れを書いていきます。

を決める

まずとして有効な値を決め、初期値とします。

左図の赤線をこのとします。

におけるの接線

次にで接するの接線を考えます。

左図の青い線をこの接線とします。

接線と軸との交点の座標が

この接線と軸との交点における座標をとします。

におけるの接線

今度は、で接するの接線を考えます。

接線と軸との交点の座標が

先ほどと同じように、接線と軸との交点における座標を、今度はとします。

におけるの接線

以下、同様に

と、を求めていきます。

接線と軸との交点の座標が、との比較

すると、やがてはと軸との交点に近づいていきます。

左図でとを比較すると、解であると軸との交点に近づいていることがわかります。

これは、がの解に収束することを意味します。

以上のようにして、ニュートン・ラフソン法での解を近似値として求めることができます。

逆数を求める

ニュートン・ラフソン法で逆数を求める際の式、

逆数を求めたい数を、Sourceよりとすることとします。を満たすときのの値が解となるようにするので、

（は逆数を求めたい数）

となります。

におけるの接線

まず接線の傾きは、

となります。

接線のy切片の座標をとすると、における接線を表す式（tangent lineよりとします）は、

（はにおける接線）
（は接線のy切片の座標）

となります。

との交点

を求めます。接点座標は、より、となります。これをに当てはめると、

（は接線のy切片の座標）
（は逆数を求めたい数）

と、が求まりました。

以上より、最終的ににおける接線は、

（はにおける接線）
（は逆数を求めたい数）

となります。

からを求める漸化式

からを求めるには、における接線のx切片の座標を求めます。

の解が

要するに、先程求めた接線がになる解を求めればいいのです。

これが漸化式となります。

あとは、この漸化式を使ってからへの変化が極僅かとなるまで収束させ続けます。

ニュートン・ラフソン法で逆数を求めるプログラムを書いてみる

Golang（Go言語、正式には単に「Go」らしいです）の入門も兼ねてGoで書いてみました。

A Tour of Go（日本語版）とPackages - The Go Programming Languageに入門者向けの読み物と、パッケージ（API？）一覧がありますので、これを読みながら書きました。

をどうするか

ニュートン・ラフソン法では、が解に近いほど早く収束します。故にはなるべく解に近づけます。

逆数を求める場合は、（逆数を求めたい数）がのときが解より大きすぎると、となってしまい、収束しないはずです。故にを満たすようにします。（のときは不等号を逆にして考えてください。）

最善策がよくわからないので、まずとし、が未満になるまでにをかけ続けるようにしました。

終了条件はどうするか

がに近づき続ける間は、続けるようにしました。

ニュートン・ラフソン法で逆数を求める関数

func reciprocal(Src float64) float64 {
	if Src == 0.0 {
		return math.NaN()
	}

	var newdist /*終了判定用変数*/, newX /*Xn+1*/ float64
	var olddist /*終了判定用変数*/, oldX /*Xn*/ float64 = math.Inf(0), Src * 0.1

	//X0を決定。Src*X0が1より小さくなるまでX0を0.1倍し続ける。
	for {
		if Src*oldX <= 1.0 {
			break
		}
		oldX *= 0.1
	}

	for {
		//漸化式でXn+1を求める
		newX = oldX * (2.0 - oldX*Src)

		//前回よりもSrc*Xn+1が1に近づいたら続行
		newdist = math.Abs(1.0 - newX*Src)
		if newdist >= olddist {
			break
		}
		olddist = newdist

		//Xnを更新
		oldX = newX
	}

	return oldX
}

ニュートン・ラフソン法で逆数を求め除算をするプログラム

上記関数を使って除算をするプログラムです。

package main

import (
	"fmt"
	"math"
	"os"
	"strconv"
)

func main() {
	var v1, v2 float64
	var stat bool

	fmt.Printf("入力1>> ")
	stat, v1 = readValidValue()
	if !stat {
		fmt.Printf("数値をいれよ\n")
		return
	}

	fmt.Printf("入力は%Fであるな\n", v1)

	fmt.Printf("入力2>> ")
	stat, v2 = readValidValue()
	if !stat {
		fmt.Printf("数値をいれよ\n")
		return
	}

	fmt.Printf("入力は%Fであるな\n", v2)

	fmt.Printf("%F/%Fは%Fである\n", v1, v2, divide(v1, v2))
}

//数値入力（ふざけた入力は受け付けない）
func readValidValue() (bool, float64) {
	var buf /*直接入力するバッファ*/, num /*入力を溜め込み、後に[]byte→string→float64と変換*/ []byte
	var dotAppeared, invalid = false, false
	var ret float64

	buf = make([]byte, 1)
	num = make([]byte, 0, 127)

	for {
		//一文字（1byte）入力
		os.Stdin.Read(buf)
		if buf[0] == byte('\n') {
			//改行コードなら入力終了
			if len(num) == 0 {
				//このとき、入力が改行コードのみだった場合、不正な入力とする
				invalid = true
			}
			break
		} else if buf[0] == byte('-') && len(num) == 0 {
			//負の数を入力するため、先頭のみ「-」を許容
			num = append(num, buf[0])

		} else if buf[0] == byte('-') && len(num) != 0 {
			//先頭以外で-を入力したら不正な入力とする
			invalid = true
			readForNewLine()
			break

		} else if buf[0] == byte('.') && len(num) >= 1 && !dotAppeared {
			//少数入力のため、２文字目以降の初めての「.」のみ許容する
			dotAppeared = true
			num = append(num, buf[0])

		} else if buf[0] == byte('.') && (len(num) == 0 || dotAppeared) {
			//先頭または二度目の.なら不正な入力とする
			invalid = true
			readForNewLine()
			break

		} else if buf[0] >= byte('0') && buf[0] <= byte('9') {
			//0〜9は常に許容
			num = append(num, buf[0])

		} else {
			//あとはすべて不正な入力とする
			invalid = true
			readForNewLine()
			break

		}
	}

	ret, _ = strconv.ParseFloat(string(num), 64)

	return !invalid, ret
}

//読み捨て用の関数
func readForNewLine() {
	buf := make([]byte, 1)

	for {
		os.Stdin.Read(buf)
		if buf[0] == byte('\n') {
			break
		}
	}
}

//除算を逆数との掛け算として求める関数
func divide(a, b float64) float64 {
	if b == 0.0 {
		return math.NaN()
	}

	return a * reciprocal(b)
}

//ニュートン・ラフソン法で逆数を求める関数
func reciprocal(Src float64) float64 {
	if Src == 0.0 {
		return math.NaN()
	}

	var newdist /*終了判定用変数*/, newX /*Xn+1*/ float64
	var olddist /*終了判定用変数*/, oldX /*Xn*/ float64 = math.Inf(0), Src * 0.1

	//X0を決定。Src*X0が1より小さくなるまでX0を0.1倍し続ける。
	for {
		if Src*oldX < 1.0 {
			break
		}
		oldX *= 0.1
	}

	for {
		//漸化式でXn+1を求める
		newX = oldX * (2.0 - oldX*Src)

		//前回よりもSrc*Xn+1が1に近づいたら続行
		newdist = math.Abs(1.0 - newX*Src)
		if newdist >= olddist {
			break
		}
		olddist = newdist

		//Xnを更新
		oldX = newX
	}

	return oldX
}

実行結果

[アホのPC]$go run main.go
入力1>> 8
入力は8.000000であるな
入力2>> 32
入力は32.000000であるな
8.000000/32.000000は0.250000である

無事に計算できてますが、特にGolangでなくてもいいような出来になりました。

geeksinhitachiprovince 2019-01-08 02:22 読者になる

広告を非表示にする

• もっと読む

コメントを書く