高校数学をプログラミングで解く（数学III編）「5-3 第ｎ次導関数」

割引あり

MK's papa

2024年11月20日 05:25

マガジンリスト　＞　数学Ⅲ編 5.微分法　＞　5-3 第ｎ次導関数

はじめに

今回は、数学IIIで学ぶ「第ｎ次導関数」について、第２次導関数をコンピュータで計算することを考えてみます。

第ｎ次導関数

まず、第ｎ次導関数について解説しておきます。

第ｎ次導関数

関数$${y=f(x)}$$を$${n}$$回微分して得られる関数。

第２次導関数の数値微分

今回は、上記で説明した第ｎ次導関数のうち、第２次導関数に注目します。第２次導関数を数値微分で求めたものと解析的に求めたものとをグラフを描いて比較してみます。

問題
次の関数の第２次導関数を求め、その数値微分と比較せよ。

$$
y=x^4-5x^3+1
$$

第２次導関数の数値微分の式

まず、第２次導関数を数値微分で求める方法を考えてみます。
第２次導関数は微分の定義に戻ってみると、

$$
y'' = \lim_{h \to 0} \frac{f'(x+h)-f'(x)}{h}
$$

と、第１次導関数$${f'(x)}$$で表すことができます。第１次導関数は

$$
f'(x) = \lim_{h \to 0} \frac{f(x+h)-f(x)}{h}
$$

と定義されますので、この式を代入することで第２次導関数の数値微分の式が求められそうです。
記事『高校数学をプログラミングで解く（数学II編）「5-1 微分係数、導関数」』で数値微分の方法として、前方差分、後方差分、中心差分の３つの方法を紹介しました。今回はこれら３つの方法のうち、中心差分を用いて第２次導関数の数値微分の式を表してみます。微分を中心差分に置き換えていくと、

$$
y'' \approx \frac{f'(x+\frac{h}{2})-f'(x-\frac{h}{2})}{h} \approx \frac{\frac{f(x+\frac{h}{2} + \frac{h}{2})-f(x+\frac{h}{2} - \frac{h}{2})}{h}-\frac{f(x-\frac{h}{2} + \frac{h}{2})-f(x-\frac{h}{2} - \frac{h}{2})}{h}}{h}
$$

となり、整理すると、

$$
y'' \approx \frac{f(x+h)-2f(x)+f(x-h)}{h^2}
$$

となります。

アルゴリズム設計

第２次導関数の数値微分の式を中心差分を用いて表すことができたので、あとは、解析的に求めた第２次導関数を求めてグラフを描いて比較します。そのために、問題の第２次導関数を求めておきます。

$$
y' = 4x^3-15x^2, \ \ y'' = 12 x^2 -30 x
$$

プログラム

それでは、第２次導関数の数値解と解析解を比較するためのプログラムを作成していきます。

今回は、記事『高校数学をプログラミングで解く（数学III編）「5-2 いろいろな関数の導関数」』で作成したスケッチ「check_derivative_sin」を再利用して作成していきます。以下の zip ファイルをダウンロードして展開または解凍してください。

そして、スケッチの名前（フォルダ名）を「check_second_derivative」と変更し、またスケッチ「check_second_derivative」内の「check_derivative_sin.pde」ファイルの名前を「check_second_derivative.pde」に変更します。そして、pdeファイル「check_second_derivative.pde」をダブルクリックしてスケッチ「check_second_derivative」の開発環境ウィンドウを立ち上げます。開発環境ウィンドウのタブ欄で「check_second_derivative」タブを選択し、そのテキストエリアのソースコードを以下で書き換えます。

float x_range = 5.0; // x軸の表示範囲 -x_rangeからx_rangeまで
float y_range = 20.0; // y軸の表示範囲 -y_rangeからy_rangeまで 

void setup(){
  size(500,500);
  noLoop();

  setAxes(x_range, y_range); // 座標軸の準備
  noFill();
  
  // グラフの描画範囲  
  float x_min = -x_range;
  float x_max = x_range;
  
  // 関数 f の第２次導関数（数値解）のグラフを描く
  stroke(0,255,0);  
  strokeWeight(5);
  draw_ddf_num(x_min, x_max);

  // 関数 f の第２次導関数（解析解）のグラフを描く
  stroke(255,0,0);
  strokeWeight(2);
  draw_ddf(x_min, x_max);

}

// 関数 f
float f(
  float x
){
  return pow(x,4.0)-5.0*pow(x,3.0)+1.0;
}

// 関数 f の第２次導関数（解析解）
float ddf(
  float x
){
  return 12.0*x*x-30.0*x;
}

// 関数 f の第２次導関数（数値解）
float ddf_num(
  float x
){
  float dx = 0.01;
  return (f(x+dx)-2.0*f(x)+f(x-dx))/dx/dx;
}

// 関数 f の第２次導関数（数値解）のグラフを描く関数
void draw_ddf_num(
  float x_min, // グラフの定義域の下限
  float x_max  // グラフの定義域の上限
){
  int plot_num = 2000; // グラフを描くための頂点の個数  
  
  // グラフを描画
  float x, y; // 関数の座標
  float X, Y; // キャンバス上の座標 
  beginShape();
  for(int i=1; i<plot_num; i++){
    x = x_min + (x_max - x_min) / plot_num * i; // 曲線上の点のx座標
    y = ddf_num(x); // 曲線上の点のyの値
    // キャンバス上の座標位置に換算
    X = width / 2.0 / x_range * x;
    Y = height / 2.0 / y_range * y;
    vertex(X, Y);
  }
  endShape();

}

// 関数 f の第２次導関数（解析解）のグラフを描く関数
void draw_ddf(
  float x_min, // グラフの定義域の下限
  float x_max  // グラフの定義域の上限
){
  int plot_num = 2000; // グラフを描くための頂点の個数  
  
  // グラフを描画
  float x, y; // 関数の座標
  float X, Y; // キャンバス上の座標 
  beginShape();
  for(int i=1; i<plot_num; i++){
    x = x_min + (x_max - x_min) / plot_num * i; // 曲線上の点のx座標
    y = ddf(x); // 曲線上の点のyの値
    // キャンバス上の座標位置に換算
    X = width / 2.0 / x_range * x;
    Y = height / 2.0 / y_range * y;
    vertex(X, Y);
  }
  endShape();

}

ソースコード1　第２次導関数の数値解と解析解を比較するためのプログラム

ソースコード1で、問題の式を関数 f として定義し、その関数 f の第２次導関数について解析的に求めたものを ddf 関数、数値微分により求めたものを ddf_num 関数として定義しています。なお、数値微分による第２次導関数 ddf_num 内で利用する差分$${ \Delta x}$$は変数 dx = 0.01 で表しています。そして、数値微分による第２次導関数 ddf_num と解析的に求めた第２次導関数 ddf のグラフを描くための関数をそれぞれ draw_ddf 関数と draw_ddf_num 関数として準備し、これらの関数を setup 関数内で呼び出すことでグラフを描いています。

ソースコード1を、スケッチ「check_second_derivative」の「check_second_derivative」タブのテキストエディタ部分に書いて実行すると、図1のように、実行ウィンドウのキャンバス上に数値微分による第２次導関数のグラフが緑色の太線で、解析的に求めた第２次導関数のグラフが赤色の線でそれぞれ描かれます。

考察

図1を見ると、数値微分による第２次導関数のグラフと解析的に求めた第２次導関数のグラフが重なっているように見えます。これは数値微分による関数と解析解が一致しているだろうということを示唆しています。
ただ、図1の極小値付近をよく見てみると、数値微分による第２次導関数のグラフ（緑色）がギザギザしているのがわかります。実際、今回の数値微分による第２次導関数 ddf_num の値を極小値付近（$${x = 1.25}$$）で見てみると、

x= 1.230 ddf_num= -18.7397
x= 1.235 ddf_num= -18.7397
x= 1.240 ddf_num= -18.749237
x= 1.245 ddf_num= -18.754005
x= 1.250 ddf_num= -18.749237
x= 1.255 ddf_num= -18.754005
x= 1.260 ddf_num= -18.749237
x= 1.265 ddf_num= -18.744469
x= 1.270 ddf_num= -18.734932

となっており、$${x = 1.25}$$で値が少し大きくなっていることがわかります。
また、数値微分による第２次導関数 ddf_num 内で利用する差分$${ \Delta x}$$をもう少し小さく dx = 0.001 としてグラフを描くと、図2のように、数値微分による第２次導関数のグラフが明らかにギザギザになり、数値微分がうまくいっていないことがわかります。

桁落ち

なぜこのようなギザギザしたグラフになるのでしょうか。原因は、桁落ちという現象に起因しています。

float型の浮動小数点数の有効桁数は 7 桁程度となっています。つまり、float型の値が

$$
0.123456789 \times 10^n
$$

となった場合に、おおよそ小数点第7位まではある程度信用できる値とみなせますが、小数点第8位以降はあまり信用できない値になります。

これを踏まえて、関数$${f(x) = x^4-5x^3+1}$$の$${x=1.25}$$付近での値を見てみます。ここで、$${\Delta x = 0.01}$$とします。

$${ f(x - \Delta x) = -0.6168906 \times 10^1 }$$
$${ f(x) = -0.63242188 \times 10^1 }$$
$${ f(x + \Delta x) = -0.6481406 \times 10^1 }$$

これらの結果を用いて、数値微分による第２次導関数の分子を計算してみると、

$${ f(x + \Delta x) -2 f(x) + f(x - \Delta x) = -0.00018744 \times 10^1 }$$

となります。$${f(x - \Delta x), f(x), f(x + \Delta x)}$$に対する信頼できる桁数は小数第7位まででしたので、$${ f(x + \Delta x) -2 f(x) + f(x - \Delta x) }$$についても小数第7位までが信頼できる値となりますが、小数第3位までは$${0}$$となっています。そのため、$${ f(x + \Delta x) -2 f(x) + f(x - \Delta x) }$$の値の小数第4位から小数第7位までの$${1874}$$の部分が信頼できる値となります。つまり、有効桁数は7桁から4桁に下がってしまいます。このような現象を桁落ちと呼んでいます。

先程、極小値付近（$${x = 1.25}$$）で数値微分による第２次導関数 ddf_num の値を見てみましたが、あらためて見てみると、

x= 1.240 ddf_num= -18.749237
x= 1.245 ddf_num= -18.754005
x= 1.250 ddf_num= -18.749237
x= 1.255 ddf_num= -18.754005
x= 1.260 ddf_num= -18.749237

のように、桁数が4桁から5桁目のあたりでおかしな値が出てきていることがわかります。つまり、有効桁数以下の値がおかしな挙動に影響していることが示唆されます。

今度は、$${\Delta x = 0.001}$$としてみます。

$${ f(x - \Delta x) = -0.6308602 \times 10^1 }$$
$${ f(x) = -0.63242188 \times 10^1 }$$
$${ f(x + \Delta x) = -0.63398542 \times 10^1 }$$

これらの結果を用いて、数値微分による第２次導関数の分子を計算してみると、

$${ f(x + \Delta x) -2 f(x) + f(x - \Delta x) = -0.00000186 \times 10^1 }$$

となります。今度は、小数第5位までが$${0}$$となっています。そのため、$${ f(x + \Delta x) -2 f(x) + f(x - \Delta x) }$$の値の小数第6位から小数第7位までの$${18}$$の部分が信頼できる値となります。つまり、有効桁数は7桁から2桁に下がってしまいます。
そして、極小値付近（$${x = 1.25}$$）で数値微分による第２次導関数 ddf_num の値を見てみると、

x= 1.240 ddf_num= -19.073486
x= 1.245 ddf_num= -19.073486
x= 1.250 ddf_num= -18.59665
x= 1.255 ddf_num= -19.550323
x= 1.260 ddf_num= -19.550323

のように、桁数が2桁から3桁目のあたりでおかしな値が出てきていることがわかります。この結果からも、有効桁数以下の値がおかしな挙動に影響していることが示唆されます。

改善案

今回のケースで、この桁落ちを回避する方法を考えてみます。

① 倍精度浮動小数点数 double 型を利用する
これまで浮動小数点数は float 型として扱ってきました。この float 型の浮動小数点数は32bitで表される数で単精度浮動小数点数とも呼ばれ、有効桁数としては7桁から8桁まででした。実は、浮動小数点数を扱う型として、倍精度浮動小数点数 double 型というものがあります。この double 型の浮動小数点数は64bitで表される数で、有効桁数としては15桁から16桁程度あります。そのため、今回の問題であれば、桁落ちしたとしても有効桁数は10桁程度は残りますので、グラフを描くためには十分な精度が保てます。実際、第２次導関数を計算する際に、float 型の変数を double 型に置き換えて計算してみます。

float x_range = 5.0; // x軸の表示範囲 -x_rangeからx_rangeまで
float y_range = 20.0; // y軸の表示範囲 -y_rangeからy_rangeまで 

void setup(){
  size(500,500);
  noLoop();

  setAxes(x_range, y_range); // 座標軸の準備
  noFill();
  
  // グラフの描画範囲  
  double x_min = -x_range;
  double x_max = x_range;
  
  // 関数 f の第２次導関数（数値解）のグラフを描く
  stroke(0,255,0);  
  strokeWeight(5);
  draw_ddf_num(x_min, x_max);

  // 関数 f の第２次導関数（解析解）のグラフを描く
  stroke(255,0,0);
  strokeWeight(2);
  draw_ddf(x_min, x_max);

}

// 関数 f
double f(
  double x
){
  return x*x*x*x-5.0*x*x*x+1.0;
}

// 関数 f の第２次導関数（解析解）
double ddf(
  double x
){
  return 12.0*x*x-30.0*x;
}

// 関数 f の第２次導関数（数値解）
double ddf_num(
  double x
){
  double dx = 0.001;
  return (f(x+dx)-2.0*f(x)+f(x-dx))/dx/dx;
}

// 関数 f の第２次導関数（数値解）のグラフを描く関数
void draw_ddf_num(
  double x_min, // グラフの定義域の下限
  double x_max  // グラフの定義域の上限
){
  int plot_num = 2000; // グラフを描くための頂点の個数  
  
  // グラフを描画
  double x, y; // 関数の座標
  float X, Y; // キャンバス上の座標 
  beginShape();
  for(int i=1; i<plot_num; i++){
    x = x_min + (x_max - x_min) / plot_num * i; // 曲線上の点のx座標
    y = ddf_num(x); // 曲線上の点のyの値
    // キャンバス上の座標位置に換算
    X = (float)(width / 2.0 / x_range * x);
    Y = (float)(height / 2.0 / y_range * y);
    vertex(X, Y);
  }
  endShape();

}

// 関数 f の第２次導関数（解析解）のグラフを描く関数
void draw_ddf(
  double x_min, // グラフの定義域の下限
  double x_max  // グラフの定義域の上限
){
  int plot_num = 2000; // グラフを描くための頂点の個数  
  
  // グラフを描画
  double x, y; // 関数の座標
  float X, Y; // キャンバス上の座標 
  beginShape();
  for(int i=1; i<plot_num; i++){
    x = x_min + (x_max - x_min) / plot_num * i; // 曲線上の点のx座標
    y = ddf(x); // 曲線上の点のyの値
    // キャンバス上の座標位置に換算
    X = (float)(width / 2.0 / x_range * x);
    Y = (float)(height / 2.0 / y_range * y);
    vertex(X, Y);
  }
  endShape();

}

ソースコード2　ソースコード1の float 型変数を double 型に置き換えたプログラム

スケッチ「check_second_derivative」の「check_second_derivative」タブのテキストエディタ部分をソースコード2に置き換えて実行すると、図3のように、実行ウィンドウのキャンバス上に数値微分による第２次導関数のグラフが緑色の太線で、解析的に求めた第２次導関数のグラフが赤色の線でそれぞれ描かれ、今度は緑色の線もギザギザにならないことがわかります。

② 桁落ちを起こすような計算をコンピュータにさせない
最善の策は桁落ちを起こすような計算をコンピュータにさせないことです。つまり、今回の問題では第２次導関数を解析的に求めることができますので、グラフを描くにしても解析的に求めた式を利用して描くようにすればいいです。
桁落ちは、値がほぼ等しく丸め誤差を持つ数値同士の減算を行った場合に起きます。例えば、$${ \sqrt{a+1}-\sqrt{a} }$$のような式も$${a}$$が大きいときに桁落ちの影響を受けそうです。ただこのような式も

$$
\sqrt{a+1}-\sqrt{a} = \frac{1}{ \sqrt{a+1}+\sqrt{a} }
$$

のように変形することで、桁落ちを避けることができます。

まとめ

今回は、数学IIIで学ぶ「第ｎ次導関数」について、第２次導関数をコンピュータで計算することを考えました。特に、関数の第２次導関数を数値微分を利用して計算しましたが、そこでは桁落ちという問題が起こることをみました。
桁落ちは、値がほぼ等しく丸め誤差を持つ数値同士の減算を行った場合に起きます。もし、計算結果がおかしかったり、グラフが思ったように描けなかったりする場合はこの桁落ちを疑ってみることも必要になるかもしれません。コンピュータでの計算ではこのような桁落ちが起こるということを覚えておいてください。
また、基本的には、桁落ちするような計算はコンピュータにさせないことをお勧めします。コンピュータで計算させる前に、桁落ちするような箇所は改善しておきましょう。桁落ちを避ける方法はいろいろありますので、その都度考えたり調べたりしてみてください。大概の桁落ちは避けられるはずです。

読んだ感想などをお寄せください

本記事を読んだ感想や質問などを以下のお問い合せフォームからお寄せください。感想、質問をいただいた方には本記事の演習問題の解答をプレゼントします。（お問合せフォームの本文に、本記事のタイトルを入れてください。）

参考文献

改訂版教科書傍用スタンダード・オリジナル数学III（数研出版、ISBN9784410209567）

演習問題

上記では、関数$${y=x^4-5x^3+1}$$の第２次導関数について、中心差分を用いました。今度は、前方差分や後方差分を用いた計算を行ってみてください。そして第２次導関数の解析的な結果と比較してみてください。

演習問題の解答例

ここから先は

1,616字 / 4画像

ログイン

この記事が気に入ったらチップで応援してみませんか？