中心差分法とは
微分とは何か聞いたことがありますか？
今回は微分について、知っていることを前提とした記事なので、微分を知らない方は基礎知識を身に着けてから、この記事を読むと面白いと思われます。
今回の内容はコンピュータで微分の問題を解く、といったものです。
まず人間が手計算などで求める微分は解析と呼ばれたりします。（解析的に解く）
解析的に解くことができる微分は連続であることが前提です。
次にコンピュータを使って微分をすることを数値微分と言います。
数値微分はコンピュータを使っているので、連続量を表すことができず、離散的になります。
そこで近似といった技を使うことで、解析的に解く微分と同じような値を得よう、といったものです。
微分の定義式は次の通りです。
微分の定義式
\begin{align}
y &= f(x)\tag{1}\\
y' &= \lim_{h\to0}\frac{f(x+h) - f(x)}{h}\tag{2}
\end{align}
補足
微分係数$f'(a)$を求めたい場合は、
\begin{align}
y'(a) = f'(a) = \lim_{h\to0}\frac{f(a+h) - f(a)}{h}
\end{align}
ちなみに、(2)式は(3)式のようにも表すことができ、これを変形することで(2)式を導出することができます。
一般的には(2)式を利用することが多いから、こちらを今回は利用します。
\begin{align}
y' &= \lim_{z\to x} \frac{f(z) - f(x)}{z - x}\tag{3}
\end{align}
$z = x + h$とすると、$z\to x$より$h\to0$と置き換えられます。したがって、(3)式は
\begin{align}
y' &= \lim_{h\to0}\frac{f(x+h) - f(x)}{x+h-x}\\
&= \lim_{h\to0}\frac{f(x+h) - f(x)}{h} \tag{2}
\end{align}
と書き表すことができるのです。
(2)式のように
\frac{f(x+h) - f(x)}{h}\tag{4}
で求める微分を前進差分法と言います。
この前進微分は解析的に解く微分の定義式と同じ形ですが、数値微分（コンピュータが近似的に解く微分）で解く場合には精度が劣ります。
そこで、さらに精度を上げた中心差分法というものをここでは使います。
中心微分は次のように表します。
中心差分の近似式
\begin{align}
\frac{f(x+h) - f(x-h)}{2h}\tag{5}
\end{align}
プログラムで実装する
簡単な微分
今回はC言語を使って数値微分を実装していきます。
ここでは例として、
\begin{align}
f(x) = x^2\\
f'(x) = 2x
\end{align}
としての計算を考えてプログラムを作成します。
また、ここで作成するコードでは微分係数を求めるものとなっているので、プログラム実行時に微分係数を求めたい$x = a$の値を入力するようにしています。
#include <stdio.h>
#include <math.h> // 必要な場合は使用する
double f(double x){
/* ここに求めたい関数を記述 */
return x*x; // 今回は y = x^2 とする
}
int main(void){
/* 微分係数を求めたい点 */
double a; // x = a のときの微分係数
printf("aの値を入力 : ");
scanf("%lf", &a);
/* 微小変化量 */
double h = 1e-5; // h = 1.0 * 10^(-5)
/* 中心差分の近似式 */
double result = (f(a+h) - f(a-h)) / (2.0*h);
/* 結果の表示 */
printf("f'(%lf) = %lf\n", a, result);
return 0;
}
実行するには、
$ gcc centDiff.c -lm
$ ./a.out
と入力します。
実行できると、aの値を入力 : と聞かれますので、$f(x) = x^2$を微分した$f'(x) = 2x$の$x$に代入したい値を入力します。
例えば、$a=3$の場合、解析的に解いたら$f'(3) = 2\times3 = 6$となります。
つまり、数値微分でも同じ値が出てこればいいわけです。検証してみると、
aの値を入力 : 3.0
f'(3.000000) = 6.000000
しっかり、同じ値が出てきていますね。
三角関数の微分
では、三角関数を使った関数の微分もやってみましょう！
今回は少し複雑な$f(x) = \sin{x} + \cos{x}$でやってみます。
ちなみに微分すると、$f'(x) = \cos{x} - \sin{x}$となります。
三角関数を使用するには、C言語で
#include <math.h>
を最初に記述する必要があります。
先程のcentDiff.cには既に記述済みですので、関数f(double x)についてだけ変更すればよいです。
また、度数法（$30^\circ$や$90^\circ$など）ではなく、弧度法（$\frac{\pi}{6}$や$\frac{\pi}{2}$など）で表す必要があるので、度数法から弧度法に変換するコードを記述する必要があります。
変換方法は簡単です。
度数法で表す角度を$x$、弧度法で表したい角度を$\theta$としたとき、
\begin{align}
180 : x &= \pi : \theta\tag{6}
\end{align}
という比率で考えることができます。
したがって、
\begin{align}
180 : x &= \pi : \theta\\
180\;\theta &= \pi x\\
\theta &= \frac{\pi x}{180}
\end{align}
という関係式になります。
これをC言語のプログラムでも実装すればよいわけです。
度数法を弧度法に変換する式
度数法の角度$x$と弧度法にしたい角度$\theta$があるとき、
\begin{align}
\theta &= \frac{\pi x}{180}\tag{7}
\end{align}
の関係が成り立ちます。
#include <stdio.h>
#include <math.h> // 今回は必要
double f(double theta){
/* ここに求めたい関数を記述 */
return sin(theta) + cos(theta); // 今回は y = sin(x) + cos(x) とする
}
int main(void){
/* 微分係数を求めたい点(度数法で受け取る) */
double a;
printf("aの値を入力 : ");
scanf("%lf", &a);
/* 弧度法に変換（例:60度ならばM_PI/3） */
double x = (M_PI * a) / 180.0;
/* 微小変化量 */
double h = 1e-5; // h = 1.0 * 10^(-5)
/* 中心差分の近似式 */
double result = (f(x+h) - f(x-h)) / (2.0*h);
/* 結果の表示 */
printf("f'(%lf) = %lf\n", a, result);
return 0;
}
実行するには、
$ gcc centDiffTrig.c -lm
$ ./a.out
と入力します。
入力値$a$は度数法です。試しに$90^\circ$で確かめてみましょう。
\begin{align}
f'(90^\circ) = f'(\pi/2) = \cos{\frac{\pi}{2}} - \sin{\frac{\pi}{2} = 0 - 1 = -1}
\end{align}
プログラムでも$-1$が出力されれば成功です。
実行結果は以下の通りです。
aの値を入力 : 90
f'(90.000000) = -1.000000
$60^\circ$の場合も試してみましょう。
aの値を入力 : 60
f'(60.000000) = -0.366025
しっかりと動作していることが分かります。
合成関数の微分
合成関数の微分についてわからない方は、以前こちらの記事を作成しましたので参考にしてください。
【参考】合成関数の微分法と偏微分法（高校・大学数学）/ 抗原抗体反応
最後に合成関数の微分をプログラムでできるのか確かめましょう。
まず、今回の例として
\begin{align}
f(x) &= \log(x^2+2)\tag{8}\\
f(x) &= \sin(\cos{x})\tag{9}
\end{align}
この2式を使って検証していきます。
また、$\log$は自然対数（$\ln$）です。
(8)式をグラフで示すとこんな感じです。
(9)式をグラフで示すとこんな感じです。
ちなみに(8)式と(9)式をそれぞれ微分すると、
(8)式は
\begin{align}
f(x) &= \log(x^2+2)\tag{8}\\
f'(x) &= \frac{2x}{x^2+2}
\end{align}
(9)式は
\begin{align}
f(x) &= \sin(\cos{x})\tag{9}\\
f'(x) &= -\sin{x}\;{\cos(\cos{x})}
\end{align}
という導関数を得ます。
それでは、まずは(8)式を数値微分するプログラムから作成していきます。
ベースとなるのはcentDiff.cです。
関数f(double x)の部分を書き換えるだけで完成するので、非常にシンプルで簡単です。
#include <stdio.h>
#include <math.h> // 今回は使用する
double f(double x){
/* ここに求めたい関数を記述 */
return log(x*x + 2); // 今回は y = log(x^2+2) とする
}
int main(void){
/* 微分係数を求めたい点 */
double a; // x = a のときの微分係数
printf("aの値を入力 : ");
scanf("%lf", &a);
/* 微小変化量 */
double h = 1e-5; // h = 1.0 * 10^(-5)
/* 中心差分の近似式 */
double result = (f(a+h) - f(a-h)) / (2*h);
/* 結果の表示 */
printf("f'(%lf) = %lf\n", a, result);
return 0;
}
コンパイル方法は下記の通りです。
$ gcc centDiff_ComFunc_1.c -lm
$ ./a.out
まずは、例として$a=2$のときの微分係数を求めることにします。
普通に手計算すると、
\begin{align}
f(x) &= \log(x^2+2)\\
f'(x) &= \frac{2x}{x^2+2}\\
f'(2) &= \frac{2\times2}{2^2+2}\\
&\simeq 0.666667
\end{align}
このような結果になります。
今度は作成したプログラムで実行してみましょう。
aの値を入力 : 2
f'(2.000000) = 0.666667
凄く正確ですね！
では、次は(9)式の場合で試してみましょう。
先に数値微分するプログラムから作成していきます。
#include <stdio.h>
#include <math.h> // 今回は使用する
double f(double theta){
/* ここに求めたい関数を記述 */
return sin(cos(theta)); // 今回は y = sin(cos(x)) とする
}
int main(void){
/* 微分係数を求めたい点 */
double a; // x = a のときの微分係数
printf("aの値を入力 : ");
scanf("%lf", &a);
/* 弧度法に変換 */
double x = (M_PI * a) / 180.0;
/* 微小変化量 */
double h = 1e-5; // h = 1.0 * 10^(-5)
/* 中心差分の近似式 */
double result = (f(x+h) - f(x-h)) / (2*h);
/* 結果の表示 */
printf("f'(%lf) = %lf\n", a, result);
return 0;
}
今回は例として$a=90^\circ$のときの微分係数を求めることにします。
普通に手計算すると、
\begin{align}
f(x) &= \sin(\cos{x})\\
f'(x) &= -\sin{x}\;{\cos(\cos{x})}\\
f'(\frac{\pi}{2}) &= -\sin{\frac{\pi}{2}}\cdot \cos(\cos{\frac{\pi}{2}})\\
&= -1\cdot\cos{0}\\
&= -1\cdot 1\\
&= -1
\end{align}
このような計算結果になります。
今度は作成したプログラムで実行してみましょう。
コンパイル方法は下記の通りです。
$ gcc centDiff_ComFunc_1.c -lm
$ ./a.out
では、実行してみましょう。
aの値を入力 : 90
f'(90.000000) = -1.000000
しっかりと動いていますね。
合成関数は手計算するより楽かも...？しれませんね。
まとめ
C言語で微分（微分係数を求める）を実装するには数値微分を使用します。
特に今回のような中心差分法などの精度を重視したアルゴリズムを使用することで結果はより、解析的に解いたときと同じ値に近似します。
数学的な関数を使用したいときはmath.hを読み込んでプログラムを作成するとよいです。