1次元版、ニュートン法(関数の勾配が0である解を求める)

func1=function(w){
	return(w^2-2)
}

dffunc1=function(w){
	return(2*w)
}
ddffunc1=function(w){
	return(2)
}

newton102=function(w,times){
	cnt=1
	for(i in 1:times){
		w_before=w
		
		w=w-dffunc1(w)/ddffunc1(w)
		printf("cnt:%d;w=%.8f,error=%.8f\n",cnt,w,abs(dffunc1(w)))
		cnt=cnt+1
		if(abs(dffunc1(w))<=1.0e-8){
			break
		}
	}

}

> w=1
> newton102(w,1e+3)
cnt:1;w=0.00000000,error=0.00000000

> w=-1
> newton102(w,1e+3)
cnt:1;w=0.00000000,error=0.00000000

これは、間違いです。

func1=function(w){
	return(w^4)
}

dffunc1=function(w){
	return(4*w^3)
}

ddffunc1=function(w){
	return(12*w^2)
}

> w=1
> newton102(w,1e+3)
cnt:1;w=0.66666667,error=1.18518519
cnt:2;w=0.44444444,error=0.35116598
cnt:3;w=0.29629630,error=0.10404918
cnt:4;w=0.19753086,error=0.03082939
cnt:5;w=0.13168724,error=0.00913463
cnt:6;w=0.08779150,error=0.00270656
cnt:7;w=0.05852766,error=0.00080194
cnt:8;w=0.03901844,error=0.00023761
cnt:9;w=0.02601229,error=0.00007040
cnt:10;w=0.01734153,error=0.00002086
cnt:11;w=0.01156102,error=0.00000618
cnt:12;w=0.00770735,error=0.00000183
cnt:13;w=0.00513823,error=0.00000054
cnt:14;w=0.00342549,error=0.00000016
cnt:15;w=0.00228366,error=0.00000005
cnt:16;w=0.00152244,error=0.00000001
cnt:17;w=0.00101496,error=0.00000000

0に近い値になっています。最急降下法に比べずっと速いですね。

> w=-1
> newton102(w,1e+3)
cnt:1;w=-0.66666667,error=1.18518519
cnt:2;w=-0.44444444,error=0.35116598
cnt:3;w=-0.29629630,error=0.10404918
cnt:4;w=-0.19753086,error=0.03082939
cnt:5;w=-0.13168724,error=0.00913463
cnt:6;w=-0.08779150,error=0.00270656
cnt:7;w=-0.05852766,error=0.00080194
cnt:8;w=-0.03901844,error=0.00023761
cnt:9;w=-0.02601229,error=0.00007040
cnt:10;w=-0.01734153,error=0.00002086
cnt:11;w=-0.01156102,error=0.00000618
cnt:12;w=-0.00770735,error=0.00000183
cnt:13;w=-0.00513823,error=0.00000054
cnt:14;w=-0.00342549,error=0.00000016
cnt:15;w=-0.00228366,error=0.00000005
cnt:16;w=-0.00152244,error=0.00000001
cnt:17;w=-0.00101496,error=0.00000000

w=1の場合と対称的です。

func1=function(w){
	return((w-3)^2*(w+3)^2)
}

dffunc1=function(w){
	return(4*(-9*w+w^3))
}

ddffunc1=function(w){
	return(12*(-3+w^2))
}

極値はw=-3,3です。

> w=1
> newton102(w,1e+3)
cnt:1;w=-0.33333333,error=11.85185185
cnt:2;w=0.00854701,error=0.30768981
cnt:3;w=-0.00000014,error=0.00000500
cnt:4;w=0.00000000,error=0.00000000

> w=2
> newton102(w,1e+3)
cnt:1;w=5.33333333,error=414.81481481
cnt:2;w=3.97476953,error=108.09453801
cnt:3;w=3.27096267,error=22.23203709
cnt:4;w=3.03033111,error=2.21707048
cnt:5;w=3.00044939,error=0.03236329
cnt:6;w=3.00000010,error=0.00000727
cnt:7;w=3.00000000,error=0.00000000

w=1では、ダメでしたが、w=2では3に収束しました。

> w=-1
> newton102(w,1e+3)
cnt:1;w=0.33333333,error=11.85185185
cnt:2;w=-0.00854701,error=0.30768981
cnt:3;w=0.00000014,error=0.00000500
cnt:4;w=-0.00000000,error=0.00000000
> w=-2
> newton102(w,1e+3)
cnt:1;w=-5.33333333,error=414.81481481
cnt:2;w=-3.97476953,error=108.09453801
cnt:3;w=-3.27096267,error=22.23203709
cnt:4;w=-3.03033111,error=2.21707048
cnt:5;w=-3.00044939,error=0.03236329
cnt:6;w=-3.00000010,error=0.00000727
cnt:7;w=-3.00000000,error=0.00000000

w=-2では-3です。