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社で自分に合うスクールを見つけましょう．後悔はさせません！
数値積分
数値積分とは，数値解析の分野の一つで，与えられる関数の定積分の値を解析的にではなく数値的に求める求積法のことです．
数値積分は，数学的な積分（アナログ）の解析手法をコンピュータ（デジタル）の数値解析に変換する手法とも言えます．
また，数学では無理数のπの真の値を表現できますが，コンピュータでは精度の限界により無理数を正確に表現できず，誤差が発生してしまいます．
なので，その誤差をいかに小さくするのかが数値解析の重要なポイントです．
本記事では，数値積分の手法として，区分求積法とその応用を説明していきます．
区分求積法とリーマン和
区分求積法は，リーマン和（Riemann Sum）による数値積分の極限として定積分を求める方法です．
また，区分求積法は，リーマン積分（Riemann integral）と似ています．
区分求積法とリーマン積分の違いは，区分求積法は区分を有理数のみで分割可能ですが，リーマン積分は区分を有理数と無理数で分割可能なことです．
リーマン和の種類には以下があります．
- 左側リーマン和（Left Riemann Sum）
- 右側リーマン和（Right Riemann Sum）
- 中点則（Midpoint Rule）
リーマン和で数値積分するコードは以下になります．
9～12行目のf関数では，下式を計算します．
$$\int_0^1 \frac{4}{1 + x^2}dx = \pi$$
区間[0, 1]で分割する区分（この場合は長方形）の個数を1，2，4，8，…，8192と変化させた場合に，円周率\(\pi\)の真の値にどのように近づくのかを数値解析します．
|
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 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #define N 8192 double f(double x) { return 4 / (1 + x * x); } int main(void) { size_t i, n; double h, left_riemann, right_riemann, midpoint; double a = 0.0, b = 1.0; h = b - a; printf("Left Riemann Sum Right Riemann Sum Midpoint Rule\n"); for (n = 1; n <= N; n *= 2) { left_riemann = right_riemann = midpoint = 0.0; for (i = 1; i <= n; i++) { /* Left riemann sum. */ left_riemann += f(a + h * (i - 1.0)); /* Right riemann sum. */ right_riemann += f(a + h * (i)); /* Midpoint rule. */ midpoint += f(a + h * (i - 0.5)); } left_riemann *= h; right_riemann *= h; midpoint *= h; printf("%.15lf %.15lf %.15lf\n", left_riemann, right_riemann, midpoint); h /= 2.0; } return 0; } |
実行結果は以下になります．
円周率\(\pi\)の真の値（3.141592653589793...）と比較すると，左側リーマン和と右側リーマン和は誤差が大きいですが，中点則は誤差が小さいことがわかります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
$ gcc riemann.c $ a.out Left Riemann Sum Right Riemann Sum Midpoint Rule 4.000000000000000 2.000000000000000 3.200000000000000 3.600000000000000 2.600000000000000 3.162352941176470 3.381176470588235 2.881176470588235 3.146800518393943 3.263988494491089 3.013988494491089 3.142894729591689 3.203441612041389 3.078441612041389 3.141918174308560 3.172679893174975 3.110179893174975 3.141674033796337 3.157176963485654 3.125926963485654 3.141612998641847 3.149394981063752 3.133769981063752 3.141597739852815 3.145496360458282 3.137683860458282 3.141593925155547 3.143545142806921 3.139638892806921 3.141592971481231 3.142569057144072 3.140615932144073 3.141592733062653 3.142080895103359 3.141104332603359 3.141592673458007 3.141836784280684 3.141348503030684 3.141592658556845 3.141714721418772 3.141470580793773 3.141592654831557 |
区分求積法の応用：台形公式，シンプソンの公式
区分求積法の応用として台形公式 （次数1のニュートン・コーツの公式） とシンプソンの公式（次数2の閉じたニュートン・コーツの公式）を紹介していきます．
ニュートン・コーツの公式の次数が高いほど誤差が小さくなります．
台形公式（Trapezoidal Rule）とシンプソンの公式（Simpson's Rule）で数値積分するコードは以下になります．
区分求積法と同じ数式を利用しています．
|
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 |
/* * Author: Hiroyuki Chishiro * License: 2-Clause BSD */ #include <stdio.h> #define N 8192 double f(double x) { return 4 / (1 + x * x); } int main(void) { size_t i, n; double h, trapezoidal, midpoint, simpson; double a = 0.0, b = 1.0; h = b - a; trapezoidal = h * (f(a) + f(b)) / 2; printf("Trapezoidal Rule Simpson's Rule\n"); for (n = 1; n <= N; n *= 2) { midpoint = 0.0; for (i = 1; i <= n; i++) { midpoint += f(a + h * (i - 0.5)); } midpoint *= h; simpson = (trapezoidal + 2 * midpoint) / 3; printf("%.15lf %.15lf\n", trapezoidal, simpson); h /= 2; trapezoidal = (trapezoidal + midpoint) / 2; } return 0; } |
実行結果は以下になります．
台形公式よりシンプソンの公式の方が次数が高いので，円周率\(\pi\)の真の値（3.141592653589793...）に近いことがわかります．
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
$ gcc newton.c $ a.out Trapezoidal Rule Simpson's Rule 3.000000000000000 3.133333333333333 3.100000000000000 3.141568627450980 3.131176470588235 3.141592502458707 3.138988494491089 3.141592651224822 3.140941612041389 3.141592653552836 3.141429893174974 3.141592653589216 3.141551963485656 3.141592653589784 3.141582481063752 3.141592653589794 3.141590110458283 3.141592653589792 3.141592017806915 3.141592653589792 3.141592494644073 3.141592653589793 3.141592613853363 3.141592653589792 3.141592643655685 3.141592653589791 3.141592651106265 3.141592653589793 |
まとめ
C言語の数値積分を紹介しました．
具体的には，区分求積法，リーマン和，リーマン積分，台形公式，シンプソンの公式を解説しました．
数値積分を使いこなすと，数学とコンピュータの両方の面白さがわかります！
数値微分を知りたいあなたはこちらからどうぞ．
数値解析した結果をグラフで表示する時は，gnuplotがおすすめです．
C言語を独学で習得することは難しいです．
私にC言語の無料相談をしたいあなたは，公式LINE「ChishiroのC言語」の友だち追加をお願い致します．
私のキャパシティもあり，一定数に達したら終了しますので，今すぐ追加しましょう！
独学が難しいあなたは，元東大教員がおすすめするC言語を学べるオンラインプログラミングスクール5社で自分に合うスクールを見つけましょう．後悔はさせません！
