これまで,数値微分の方法を紹介してきました。今回は,そのまとめとして実際にコンピュータに計算させましょう。
問題 関数を数値微分するプログラムを作りましょう。
前進差分と中心差分それぞれの場合で求めましょう。その際,導関数から求めた値と並べて出力させてください。

解説
問題 関数を数値微分するプログラムを作りましょう。
先ずは,解析的に求めた導関数の式を示します。

(1)はともかく,(2)はいやですよね。こんなのは,導出が正しいかどうか不安でしようがありません。数値計算で求めた値が使えるなら,そうさせていただきたいところです。
それでは早速,前進差分と中心差分を適用してみましょう。微少な量hは前回検討した結果を用いて,前進差分では2-26,中心差分では2-52を用いて実行してみましょう。
ソースコード:Forward CentralDifference01.java
01: class Forward_CentralDifference01 {
02:
03: public static double f1(double x) {
04: double a = 2;
05: double b = 4;
06: return a*Math.pow(x,2)+b;
07: }
08:
09: public static double f2(double u) {
10: return Math.sqrt(1/(1+Math.pow(Math.E,-2*u)));
11: }
12:
13: public static void main(String[] args) {
14: double x = 2;
15: double hf = Math.pow(2,-26);
16: double hc = Math.pow(2,-52);
17:
18: double fdiff1 = 0; //forward difference 前進差分
19: double cdiff1 = 0; //central difference 中心差分
20: System.out.println("f1(x) = " + f1(x));
21: System.out.println("f1(x+hf) = " + f1(x+hf));
22: fdiff1 = ( f1(x+hf) - f1(x) )/hf;
23: System.out.println("fdiff1 = " + fdiff1);
24: cdiff1 = ( f1(x+hc) - f1(x-hc) )/(2*hc);
25: System.out.println("f1(x+hc) = " + f1(x+hc));
26: System.out.println("f1(x-hc) = " + f1(x-hc));
27: System.out.println("cdiff1 = " + cdiff1);
28: System.out.println("TrueDiff = " + 4*x);
29:
30: double fdiff2 = 0;
31: double cdiff2 = 0;
32: System.out.println("f2(x) = " + f2(x));
33: System.out.println("f2(x+hf) = " + f2(x+hf));
34: fdiff2 = ( f2(x+hf) - f2(x) )/hf;
35: System.out.println("fdiff2 = " + fdiff2);
36: cdiff2 = ( f2(x+hc) - f2(x-hc) )/(2*hc);
37: System.out.println("f2(x+hc) = " + f2(x+hc));
38: System.out.println("f2(x-hc) = " + f2(x-hc));
39: System.out.println("cdiff2 = " + cdiff2);
40: System.out.println("TrueDiff = "
41: + Math.pow(Math.E,-2*x)
42: / (Math.sqrt(Math.pow(1+Math.pow(Math.E,-2*x),3))) );
43:
44: }// end of main
45: }// end of class
以下はその実行結果です。
出力画面
C:>java ForwardDifference01
f1(x) = 12.0
f1(x+hf) = 12.00000011920929
fdiff1 = 8.0
f1(x+hc) = 12.0
f1(x-hc) = 11.999999999999998
cdiff1 = 4.0
TrueDiff = 8.0
f2(x) = 0.9909660892472095
f2(x+hf) = 0.9909660895128037
fdiff2 = 0.017823725938796997
f2(x+hc) = 0.9909660892472095
f2(x-hc) = 0.9909660892472095
cdiff2 = 0.0
TrueDiff = 0.01782372414651308
おやおや,(1)の中心差分の値cdiff1がおかしいですね。(2)の中心差分の値も少々疑問です。TrueDiffに示されているように,0ではないはず。微少な量 h が小さすぎて,計算途中で発生した桁落ちが大きく響いているのでしょう。
今回は微分値の真の値が分かりますから,様々な h を用いて数値計算し,真の値と比較してみましょう。最も真の値に近くなる h はどんな量なのでしょうか。
ソースコード:ForwardDifference02.java
01: class ForwardDifference02 {
02:
03: public static double f1(double x) {
04: double a = 2;
05: double b = 4;
06: return a*Math.pow(x,2)+b;
07: }
08:
09: public static double f2(double u) {
10: return Math.sqrt(1/(1+Math.pow(Math.E,-2*u)));
11: }
12:
13: public static void main(String[] args) {
14: double x = 2;
15:
16: for (int i=1 ; i<52 ; i++){
17: double h = Math.pow(2,-i);
18: double diff = (f1(x+h) - f1(x))/h;
19: double diffc = (f1(x+h) - f1(x-h))/(2*h);
20: System.out.println("i("+i+") : h("+h+") : diff = "
21: +diff+" : diffc = "+diffc);
22: }
23: for (int i=1 ; i<52 ; i++){
24: double h = Math.pow(2,-i);
25: double diff = (f2(x+h) - f2(x))/h;
26: double diffc = (f2(x+h) - f2(x-h))/(2*h);
27: System.out.println("i("+i+") : h("+h+") : diff = "
28: +diff+" : diffc = "+diffc);
29: }
30:
31: }// end of main
32: }// end of class
以下はその実行結果です。不要な部分は省略します。先ずは(1)について。
出力画面
(略)
i(25) : h(2.9802322387695312E-8) : diff = 8.000000059604645 : diffc = 8.0
i(26) : h(1.4901161193847656E-8) : diff = 8.0 : diffc = 8.0
(略)
diffは前進差分を用いたときの数値微分値です。 i が25,これは微少量 h が2-25のときまでは誤差部分を追い込んでいくことが出来ていますが, i が26になった時点で追い切れなくなっています。計算途中で桁落ちしてしまったのでしょう。倍精度実数型の精度としては半分程度のところです。diffcは中心差分を用いたときの数値微分値です。最初から最後まで,8.0を出力し続けます。今回の数値微分を行った目的位置 x =2 では,倍精度実数型で表現可能な大きさの誤差を発生しなかった,ということでしょう。

