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

第71回 微分・積分の数学 差分法 [後編]

この記事を読むのに必要な時間:およそ 2 分

これまで,数値微分の方法を紹介してきました。今回は,そのまとめとして実際にコンピュータに計算させましょう。

問題 関数を数値微分するプログラムを作りましょう。

前進差分と中心差分それぞれの場合で求めましょう。その際,導関数から求めた値と並べて出力させてください。

解説

問題 関数を数値微分するプログラムを作りましょう。

先ずは,解析的に求めた導関数の式を示します。

(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 では,倍精度実数型で表現可能な大きさの誤差を発生しなかった,ということでしょう。

著者プロフィール

平田敦(ひらたあつし)

地方都市の公立工業高等学校教諭。趣味はプログラミングと日本の端っこ踏破旅行。2010年のLotYはRuby。結城浩氏のような仕事をしたいと妄想する30代後半♂。

コメント

コメントの記入