C言語の数値微分を教えて！
こういった悩みにお答えします．
本記事の信頼性
- リアルタイムシステムの研究歴12年．
- 東大教員の時に，英語でOS（Linuxカーネル）の授業．
- 2012年9月～2013年8月にアメリカのノースカロライナ大学チャペルヒル校（UNC）コンピュータサイエンス学部で客員研究員として勤務．C言語でリアルタイムLinuxの研究開発．
- プログラミング歴15年以上，習得している言語: C/C++，Python，Solidity/Vyper，Java，Ruby，Go，Rust，D，HTML/CSS/JS/PHP，MATLAB，Assembler (x64，ARM)．
- 東大教員の時に，C++言語で開発した「LLVMコンパイラの拡張」，C言語で開発した独自のリアルタイムOS「Mcube Kernel」をGitHubにオープンソースとして公開．
- 2020年1月～現在はアメリカのノースカロライナ州チャペルヒルにあるGuarantee Happiness LLCのCTOとしてECサイト開発やWeb/SNSマーケティングの業務．2022年6月～現在はアメリカのノースカロライナ州チャペルヒルにあるJapanese Tar Heel, Inc.のCEO兼CTO．
- 最近は自然言語処理AIとイーサリアムに関する有益な情報発信に従事．
- （AI全般を含む）自然言語処理AIの論文の日本語訳や，AIチャットボット（ChatGPT，Auto-GPT，Gemini（旧Bard）など）の記事を50本以上執筆．アメリカのサンフランシスコ（広義のシリコンバレー）の会社でプロンプトエンジニア・マネージャー・Quality Assurance（QA）の業務委託の経験あり．
- （スマートコントラクトのプログラミングを含む）イーサリアムや仮想通貨全般の記事を200本以上執筆．イギリスのロンドンの会社で仮想通貨の英語の記事を日本語に翻訳する業務委託の経験あり．
こういった私から学べます．
C言語を独学で習得することは難しいです．
私にC言語の無料相談をしたいあなたは，公式LINE「ChishiroのC言語」の友だち追加をお願い致します．
私のキャパシティもあり，一定数に達したら終了しますので，今すぐ追加しましょう！
独学が難しいあなたは，元東大教員がおすすめするC言語を学べるオンラインプログラミングスクール5社で自分に合うスクールを見つけましょう．後悔はさせません！
数値微分
数値微分とは，数値解析の分野の一つで，与えられる関数の常微分の値を解析的にではなく数値的に求める手法のことです．
数値微分は，数学的な微分（アナログ）の解析手法をコンピュータ（デジタル）の数値解析に変換する手法とも言えます．
数値微分を利用することで，微分を解析的に求めることが困難な場合に，その近似値を数値的に求めることが可能になります．
数値微分のニュートン法と二分法を知りたいあなたはこちらからどうぞ．
C言語で数値微分：前進差分，後退差分，中心差分
C言語で数値微分の前進差分，後退差分，中心差分の解法を紹介していきます．
前進差分は，一つ前の値との差分を利用して下式になります．
$$f'(x) = \frac{f(x + h) - f(x)}{h}$$
後退差分は下式になります．
$$f'(x) = \frac{f(x) - f(x - h)}{h}$$
中心差分は下式になります．
$$f'(x) = \frac{f(x + h) - f(x - h)}{2h}$$
C言語で解く数式は下式になります．
\begin{eqnarray}
f(x) &=& x^2 \\
f'(x) &=& 2x
\end{eqnarray}
前進差分，後退差分，中心差分を利用するコードは以下になります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #define H 0.0001 typedef double (*func_t)(double); double f(double x) { return x * x; } double front_num_diff(func_t f, double x) { return ((*f)(x + H) - f(x)) / H; } double back_num_diff(func_t f, double x) { return ((*f)(x) - f(x - H)) / H; } double midpoint_num_diff(func_t f, double x) { return ((*f)(x + H) - f(x - H)) / (2 * H); } int main(void) { printf("front: %lf\n", front_num_diff(f, 3)); printf("back: %lf\n", back_num_diff(f, 3)); printf("midpoint: %lf\n", midpoint_num_diff(f, 3)); return 0; } |
実行結果は以下になります．
\(f'(x)\)で\(x = 3\)の時の真の値\(f'(3) = 2 * 3 = 6\)と比較すると，前進差分は6.000100，後退差分は5.999900と誤差が発生していますが，中心差分は6.000000と誤差がほとんどないことがわかります．
※Hの値（0.0001）を小さくすると，その誤差が小さくなります．
|
1 2 3 4 5 |
$ gcc num_diff.c $ a.out front: 6.000100 back: 5.999900 midpoint: 6.000000 |
常微分方程式
常微分方程式とは，数学において未知関数とその導関数からなる等式で定義される方程式である微分方程式の一種で，未知関数が一つの変数を持つもののことを言います．
本記事では，常微分方程式の解法として，オイラー法，3次テイラー展開，4次ルンゲクッタ法を紹介していきます．
オイラー法と4次ルンゲクッタ法は，keisanの常微分方程式で解答を計算できますので，参考にして下さい．
各々のリンクは以下になります．
C言語で常微分方程式：オイラー法，3次テイラー展開，4次ルンゲクッタ法
C言語で常微分方程式のオイラー法，3次テイラー展開，4次ルンゲクッタ法を解法を紹介していきます．
解く常微分方程式は下式になります．
\begin{eqnarray*}
F(x, y) &=& 1 - y^2 \\
x_0 &=& 0 \\
y_0 &=& 0 \\
x_n &=& 1
\end{eqnarray*}
オイラー法，3次テイラー展開，4次ルンゲクッタ法を利用するコードは以下になります．
上式の常微分方程式の解答の\(y = f(x) = \tanh x\)と比較します．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #include <math.h> double square(double x) { return x * x; } double F(double x, double y) { return 1 - square(y); } double Fx(double x, double y) { return -2 * y * F(x, y); } double Fxx(double x, double y) { return -2 * (square(F(x, y)) + y * Fx(x, y)); } /* Euler Method */ double euler(double x0, double y0, double xn, int n, int nprint) { int i; double x, y, h; x = x0; y = y0; h = (xn - x0) / n; for (i = 1; i <= n; i++) { y += F(x, y) * h; x = x0 + i * h; if (i % nprint == 0) { printf("%lf %.15lf\n", x, y); } } return y; } /* 3rd-Order Taylor Method */ double taylor3(double x0, double y0, double xn, int n, int nprint) { int i; double x, y, h; x = x0; y = y0; h = (xn - x0) / n; for (i = 1; i <= n; i++) { y += h * (F(x, y) + (h / 2) * (Fx(x, y) + (h / 3) * Fxx(x, y))); x = x0 + i * h; if (i % nprint == 0) { printf("%lf %.15lf\n", x, y); } } return y; } /* 4th-Order Runge-Kutta Method */ double runge4(double x0, double y0, double xn, int n, int nprint) { int i; double x, y, h, h2, f1, f2, f3, f4; x = x0; y = y0; h = (xn - x0) / n; h2 = h / 2; for (i = 1; i <= n; i++) { f1 = h * F(x, y); f2 = h * F(x + h2, y + f1 / 2); f3 = h * F(x + h2, y + f2 / 2); f4 = h * F(x + h, y + f3); x = x0 + i * h; y += (f1 + 2 * f2 + 2 * f3 + f4) / 6; if (i % nprint == 0) { printf("%lf %.15lf\n", x, y); } } return y; } int main(void) { int i, n; int divs[] = {10, 20, 50, 100, 200, 500}; int nr_divs = sizeof(divs) / sizeof(divs[0]); for (i = 0; i < nr_divs; i++) { n = divs[i]; printf("Euler: n = %d\n", n); euler(0, 0, 1, n, n / 5); printf("3rd-Order Taylor: n = %d\n", n); taylor3(0, 0, 1, n, n / 5); printf("4th-Order Runge-Kutta: n = %d\n", n); runge4(0, 0, 1, n, n / 5); } printf("\nAnswer:\n"); for (i = 1; i <= 5; i++) { printf("%lf %.15lf\n", i / 5.0, tanh(i / 5.0)); } return 0; } |
実行結果は以下になります．
オイラー法が一番誤差が大きく，3次テイラー展開が中間，4次ルンゲクッタ法が一番誤差が小さいことがわかります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 |
$ gcc ord_diff.c -lm $ a.out Euler: n = 10 0.200000 0.199000000000000 0.400000 0.386335045740799 0.600000 0.549186870811072 0.800000 0.680706899220737 1.000000 0.780440676845054 3rd-Order Taylor: n = 10 0.200000 0.197366368474459 0.400000 0.379913111716298 0.600000 0.536988222592680 0.800000 0.663965380725370 1.000000 0.761529025004616 4th-Order Runge-Kutta: n = 10 0.200000 0.197375143874743 0.400000 0.379948535990683 0.600000 0.537048794343463 0.800000 0.664035622262367 1.000000 0.761592708599983 Euler: n = 20 0.200000 0.198260586027217 0.400000 0.383217258968393 0.600000 0.543113359463835 0.800000 0.672266095167161 1.000000 0.770854654364196 3rd-Order Taylor: n = 20 0.200000 0.197373985283189 0.400000 0.379944245133548 0.600000 0.537041838766193 0.800000 0.664027995938668 1.000000 0.761586279001855 4th-Order Runge-Kutta: n = 20 0.200000 0.197375309031945 0.400000 0.379948935322221 0.600000 0.537049518977630 0.800000 0.664036700145941 1.000000 0.761594068777302 Euler: n = 50 0.200000 0.197746258688711 0.400000 0.381272676063956 0.600000 0.539472787385898 0.800000 0.667303845935648 1.000000 0.765261460461343 3rd-Order Taylor: n = 50 0.200000 0.197375226706909 0.400000 0.379948651691020 0.600000 0.537049070266235 0.800000 0.664036214165978 1.000000 0.761593661317164 4th-Order Runge-Kutta: n = 50 0.200000 0.197375319936038 0.400000 0.379948961562139 0.600000 0.537049565773961 0.800000 0.664036768498093 1.000000 0.761594153773258 Euler: n = 100 0.200000 0.197563539088527 0.400000 0.380613422041230 0.600000 0.538260713247089 0.800000 0.665666260270956 1.000000 0.763421826145240 3rd-Order Taylor: n = 100 0.200000 0.197375308201817 0.400000 0.379948923079965 0.600000 0.537049504820536 0.800000 0.664036700976171 1.000000 0.761594094507915 4th-Order Runge-Kutta: n = 100 0.200000 0.197375320206804 0.400000 0.379948962211840 0.600000 0.537049566921649 0.800000 0.664036770157770 1.000000 0.761594155820372 Euler: n = 200 0.200000 0.197470111853204 0.400000 0.380281831122642 0.600000 0.537655017061361 0.800000 0.664850509056482 1.000000 0.762506510609003 3rd-Order Taylor: n = 200 0.200000 0.197375318701253 0.400000 0.379948957336280 0.600000 0.537049559220517 0.800000 0.664036761620056 1.000000 0.761594148298353 4th-Order Runge-Kutta: n = 200 0.200000 0.197375320223771 0.400000 0.379948962252511 0.600000 0.537049566993265 0.800000 0.664036770260985 1.000000 0.761594155947334 Euler: n = 500 0.200000 0.197413399937114 0.400000 0.380082261672787 0.600000 0.537291716608732 0.800000 0.664362025026982 1.000000 0.761958744411563 3rd-Order Taylor: n = 500 0.200000 0.197375320126594 0.400000 0.379948961939569 0.600000 0.537049566500070 0.800000 0.664036769714912 1.000000 0.761594155466589 4th-Order Runge-Kutta: n = 500 0.200000 0.197375320224875 0.400000 0.379948962255156 0.600000 0.537049566997913 0.800000 0.664036770267674 1.000000 0.761594155955550 Answer: 0.200000 0.197375320224904 0.400000 0.379948962255225 0.600000 0.537049566998035 0.800000 0.664036770267849 1.000000 0.761594155955765 |
まとめ
C言語の数値微分を紹介しました．
具体的には，前進差分，後退差分，中心差分，常微分方程式，オイラー法，3次テイラー展開，4次ルンゲクッタ法を解説しました．
数値積分を知りたいあなたはこちらからどうぞ．
数値解析した結果をグラフで表示する時は，gnuplotがおすすめです．
C言語を独学で習得することは難しいです．
私にC言語の無料相談をしたいあなたは，公式LINE「ChishiroのC言語」の友だち追加をお願い致します．
私のキャパシティもあり，一定数に達したら終了しますので，今すぐ追加しましょう！
独学が難しいあなたは，元東大教員がおすすめするC言語を学べるオンラインプログラミングスクール5社で自分に合うスクールを見つけましょう．後悔はさせません！
