【C言語】台形公式で数値積分を行う

このページではC言語における数値積分のやり方について説明していきます！このページで扱うのは、数値積分における「台形公式」を用いた手法になります。

高校や大学等の授業で積分を学んで苦手だった人も多いと思います。ですが、安心してください。まず、数値積分は、これらの授業で学んだような解析的な難しい解き方で行うものではありません。特に、今回紹介する手法においては、極論すれば台形の面積さえ求められれば積分結果を得ることが可能です。なので、特に難しい知識もなく、簡単に積分を行えるようになります。

数値積分

まずは、数値積分について説明していきます。積分や微分は奥が深く、私も詳しくないので、ここでは高校数学程度の学習内容に基づいた説明を行っていきます。また、１変数の積分を行うことを前提に解説を行っていきます。

定積分の解き方

まず、高校で学んだ定積分は、基本的に下記のような式により積分を求めていたと思います。

$$ \int_{a}^{b}f(x)dx = F(b) – F(a) $$

この、\( F(x) \) は \( f(x) \) を不定積分した結果です。また、\( a \) と \( b \) はそれぞれ積分区間の下端と上端となります。

つまり、定積分の結果は、\( f(x) \) を不定積分した結果 \( F(x) \) を求め、この \( F(x) \) における \( x = b \) の値と \( x = a \) の値の差、すなわち \( F(b) – F(a) \) を求めることで得られることになります。

で、この不定積分を求めるのが難しい…。例えば、下記の不定積分の結果は何になるでしょうか？受験勉強した方や学生さんであれば楽勝に解けるかもしれないですが、解けない方も多くおられると思います。

$$ \int \frac{1}{1 + x^2}dx $$

要は、定積分を求めるときに不定積分を求めるのが厄介です。様々な公式を駆使して求める必要があります。なんですが、数値積分の場合は不定積分を行うことなく定積分の結果を求めることが可能です。

スポンサーリンク

数値積分

では、数値積分とは何かというと、ここまで説明してきたようなやり方ではなく、もともとの関数 \( f(x) \) を数値的に評価して積分する方法のことを言います。

数値積分を理解する上で重要なポイントは、１変数の関数に対する積分は「面積を求めること」であることという点なります。より具体的には、\( f(x) \) が \( x \) の変化に伴って描く曲線（直線の場合もあり）と、積分区間 \( [a, b] \) を示す \( x = a \) と \( x = b \) の縦線と、\( x \) 軸との４つの線で囲まれる領域の面積を求めることになります。

つまり、この領域の面積は \( \int_{a}^{b}f(x)dx \) と一致します。であれば、前述のような不定積分を利用せずに、他の何らかの手段でこの領域の面積さえ求めてやれば、不定積分をしなくても \( \int_{a}^{b}f(x)dx \) の結果を求めることが出来ることになります。

ということで、面積さえ求められれば積分を求めることが出来ることになります。

簡単な例で考えてみましょう！

例えば、\( \int_{2}^{8} x dx \) の積分値を求めようと思うと、\( \int x dx = \frac{1}{2} x^2 \) と不定積分結果が得られるため、定積分は \( \frac{1}{2} 8^2 – \frac{1}{2} 2^2 = 30\) となります。

で、先ほど説明したように、この積分は面積から求めるこも可能です。\( f(x) = x \) をグラフとして描画し、さらに \( x = 2 \) と \( x = 8 \) の縦線を描画してやれば、下図の青背景部分のような台形が浮かび上がります。

そして、この台形は、上底が \( f(2) = 2 \)、下底が \( f(8) = 8 \)、高さが \( 6 \) なので、この台形の面積は下記のようにして求めることが出来ます。

$$ (2 + 8) \times 6 / 2 = 30 $$

ここで、上記の面積と、先ほど不定積分を行ってから求めた定積分の値とを比べれば、これらが一致していることが確認できると思います。このように、１変数の関数の積分値は面積から求めることが可能です。そして、面積を求める場合は、直接 \( f(x) \) の値を評価した数値 (上記であれば \( f(2) \) と \( f(8) \) )を用いて求めることから、数値積分と呼ばれます。 

ただ、先ほどは \( f(x) \) が線形に変化するため、単に台形の面積を求めるだけで積分を求めることが出来ました。ですが、\( f(x) \) が曲線を描くような場合、例えば \( f(x) = x^2 \) や、下図のように \( x \) に対して変化する \( f(x) \) に関しては面積を求めることは困難です。

なので、正確に面積を求めるのではなく、近似的に面積を求める手法をとることになります。

じゃあ、どうやって近似すれば良いでしょうか？

この近似を行う方法には下記のようなものがあります。これらは全て、面積を求めるという考え方で積分を行う手法であり、これらの違いは近似の仕方のみとなります。

• 長方形近似
• 台形公式
• シンプソンの公式

そして、今回紹介するのは、上記のうちの「台形公式」による手法となります。この台形公式においては、その名のとおり、面積を求める領域を複数の台形に近似して面積を求めることになります。

台形公式における数値積分

この台形公式における数値積分は、下記のような手順で計算を行うことになります。

面積を求める領域を分割する

まず、\( f(x) \) と \( x = a \)、\( x = b \)、さらに \( x \) 軸の４つの線で囲まれる領域に縦方向の補助線を入れ、面積を求める領域を \( x \) 軸に対して等間隔に分割します。

スポンサーリンク

分割した各領域を台形に近似する

そして、分割後の各領域の \( f(x) \) の曲線と、先ほど追加した補助線（＆積分区間の下端と上端の縦線）との交点２つを直線で結びます。つまり、各領域の中で \( f(x) \) が線形に変化すると仮定したときの線を描画します。

そうすると、\( f(x) \) と \( x = a \)、\( x = b \)、さらに \( x \) 軸の４つの線で囲まれる領域は、複数の台形で近似されることになります。近似なので誤差は含まれることになりますが、４つの線で囲まれる領域全体のおおよその面積は、分割された台形の面積の総和で求められることになります。

そして、４つの線で囲まれる領域全体（上図の右側の車線部の領域）の面積は、 \( \int_{a}^{b}f(x)dx \) の結果と一致します。したがって、先ほど近似された全ての台形の面積の総和を求めることで、\( \int_{a}^{b}f(x)dx \) のおおよその値を得ることができることになります。

各台形の面積の総和を求める

で、１つの台形の面積に関しては下記の式で求めることが出来ます。

$$ \frac{1}{2} \times (上底 + 下底) \times 高さ $$

ここで、各台形の \( x \) 軸方向に対する長さを \( h \) とすれば、１番左側に存在する台形に関しては、左側の縦線は \( x = a \) に存在し、さらに 右側の縦線は \( x = a + h \) の位置に存在することになります。さらに、左側の縦線を上底、右側の縦線を下底とすれば、上底の長さは \( f(a) \) であり、下底の長さは \( f(a + h) \) ということになりますね。

したがって、１番左側に存在する台形の面積は下記の式で求まることになります。

$$ \frac{1}{2} \times \{ f(a) + f (a + h) \} \times h $$

その右隣の台形に関しては、左側の縦線は \( x = a + h \) に存在し、さらに 右側の縦線は \( x = a + 2 \times h \) の位置に存在することになりますので、先ほどと同様にして面積が下記の式で求まることになります。

$$ \frac{1}{2} \times \{ f(a + h) + f(a + 2 \times h) \} \times h $$

一般的な形式に直せば、左から \( i \) 番目の台形の面積は下記の式によって求まることになります（１番左側の台形を \( 0 \) 番目としています）。

$$ \frac{1}{2} \times \{ f(a + i \times h) + f(a + (i + 1) \times h) \} \times h $$

さらに、分割数を \( N \) とすれば、全台形の総和は \( \sum \) を用いて下記の計算で求めることができることになります。そして、この結果は \( \int_{a}^{b}f(x)dx \) とほぼ一致することになります。

$$ \sum_{i=0}^{N-1}  \frac{1}{2} \times \{ f(a + i \times h) + f(a + (i + 1) \times h) \} \times h  $$

これは、プログラム的に言えば、\( i = 0 \) から \( i = N – 1 \) までのループの中で、前述の式で求まる台形の面積の計算を行い、その結果を足し合わせるようにしてやればよいだけです。例えばC言語であれば、下記のような for ループにより、分割した台形の面積の総和、すなわち数値積分の結果 sum を得ることが出来ることになります。

double sum = 0;
for (int i = 0; i < N; i++) {
    sum += 0.5 * (f(a + i * h) + f(a + ((i + 1) * h))) * h;
}

ということで、結局積分は台形の面積の計算で求めることができ、さらに分割数を \( N \) とすれば、\( N \) 個の台形の面積の総和により求めることが出来ます。

分割数を多くして誤差を減らす

ただし、この方法で積分を求めると誤差が発生することがあります。具体的には、\( f(x) \) が \( x \) の変化に対して線形的にではなく曲線的に変化するような場合、例えば \( f(x) \) が２次関数や３次関数のような場合に誤差が発生することになります。

この誤差は、分割した領域を台形に近似して面積を求めることが原因で発生します。本来であれば \( f(x) \) の値が \( x \) の変化に応じて曲線的に変化するのに、\( f(x) \) の値が各領域内で線形に値が変化するように無理やり各領域の形を変形させています。

そのため、\( f(x) \) が描く曲線と、台形に近似するために描いた直線の２つの線の間の領域の面積の総和分の誤差が発生することになります。

ただ、この誤差は、分割数を多くすることで小さくすることが出来ます。例えば、上の図だと誤差が大きく発生していますが、下の図のように分割数を３倍にしてやることで誤差が小さくなっていることが確認できると思います。

そして、分割数が無限大であれば誤差はゼロとなります。が、分割数が無限大ということは、面積の総和を求めるときのループ処理の回数も無限大になるため、プログラムが終了しないことになってしまいます…。なので、分割数を無限大にすることは不可能ではあるのですが、分割数を多くすることで誤差を小さくできることは覚えておいた方が良いですし、これが誤差を小さくするためのプログラム開発時のポイントとなります。

スポンサーリンク

台形公式による数値積分を行うプログラム

最後に、台形公式によるプログラムのソースコードの例をサンプルとして示していきたいと思います。

関数に対する数値積分

最初に、引数を x とする何らかの関数に対する数値積分を行うプログラムのソースコードのサンプルを示します。

#include <stdio.h>
#include <math.h>

#define N 256

typedef double (*FUNC)(double);

double f(double x) {
    return sin(x) + cos(x);
}

double integral(double a, double b, int n, FUNC f) {
    double h = (b - a) / n;

    double sum = 0;
    for (int i = 0; i < n; i++) {
        sum += 0.5 * (f(a + i * h) + f(a + ((i + 1) * h))) * h;
    }
    return sum;
}

int main(void) {
    double result = integral(0, M_PI / 2, N, f);
    printf("%lf\n", result); 
}

上記の integral が台形公式による数値積分を行う関数となります。

integral は、積分区間の下端 a と上端 b、さらに積分区間を分割する数 n、さらに積分を行う関数 f を引数で受け取り、下記の積分の結果を返却する関数となっています。

$$ \int_{a}^{b} f(x) dx $$

そして、上記のソースコードでは、各引数に下記を指定して integral 関数を実行しています。

• a : 0
• b : M_PI / 2
• n : N (256)
• f : func

M_PI は円周率 \( \pi \) を表す定義値となります。M_PI に関して詳しく知りたい方は下記ページを参照してください。

C言語で円周率πを扱う方法

また、func 関数では、sin(x) + cos(x) の計算を行っているため、上記のソースコードをコンパイルして実行すれば、下記式の計算が行われた結果が表示されることになります。そして、この計算が数値積分により求められ、その精度は分割数 n によって変化することになります。

$$ \int_{0}^{ \pi / 2} \{ \sin x  + \cos x \} dx $$

特に引数 f の型が関数ポインタとなっているのでややこしいですが、要は引数 f には「double 型の引数１つを受け取り、さらに返却値の型が double となる関数」の関数名を指定することが可能です。

この関数ポインタについては下記ページで説明していますので、関数ポインタについて詳しく知りたい方は別途下記ページを参照していただければと思います。

C言語の関数ポインタについて解説

integral 関数の説明に戻ると、この関数では、まず引数で与えられる積分区間と分割数より、面積を求める全体の領域を n 個に分割した時の各領域の幅 h を求めています。

double h = (b - a) / n;

続いて i に対する for ループを行い、このループの中で各台形の面積を求めています。

第 i 番目の台形の上底の長さ、下底の長さ、高さはそれぞれ下図のようになりますので、

第 i 番目の台形の面積は下記で求まることになります。

0.5 * (f(a + i * h) + f(a + ((i + 1) * h))) * h

そして、この各台形の面積を、ループ前に 0 で初期化した変数 sum に足し合わせていけば、ループが終わったタイミングで全台形の面積の総和が sum に格納されていることになります。この面積の総和 sum が関数 f を積分区間 [a, b] で積分した結果となりますので、その sum を返却して関数を終了しています。

double sum = 0;
for (int i = 0; i < n; i++) {
    sum += 0.5 * (f(a + i * h) + f(a + ((i + 1) * h))) * h;
}
return sum;

上記のソースコードをコンパイルして実行した場合、この sum の値は次のような値となります。

1.999994

実際に \( \int_{0}^{ \pi / 2} \{ \sin x  + \cos x \} dx \) を計算してみると下記のような式変形によって \( 2 \) が求まります。

\begin{align} \int_{0}^{ \pi / 2} \{ \sin x + \cos x \} dx &= \int_{0}^{ \pi / 2} \sin x dx + \int_{0}^{ \pi / 2} \cos x  dx \\ &= – \cos \frac{\pi}{2} + \cos 0 + \sin \frac{\pi}{2} – \sin 0 \\ &= -0 + 1 + 1 – 0 \\ &= 2 \end{align}

これらの結果より、台形の面積の総和による数値積分によって、定積分のおおよその結果が得られていることが確認できると思います。とはいえ、数値積分には誤差が含まれており、上記の場合は 0.000006 の誤差が含まれています。この誤差は、前述のとおり分割数を増やすことで小さくすることが出来ます。上記のソースコードの N の値を大きくすれば、誤差が小さくなることを確認できると思います。逆に N の値を小さくすれば誤差が大きくなることも確認できます。

また、上記のソースコードでは積分対象の関数を \( \sin x  + \cos x \) としましたが、func 関数を変更したり、新たに関数を定義してその関数を integral の引数 f に指定してやれば、もっと別の関数に対する積分を行うことも可能です。是非色々試してみていただければと思います！

配列に対する数値積分

先ほどは関数に対して数値積分を行いましたが、定期的にサンプリングした値が格納された数列に対して積分を行うようなことも可能です。

例えば、あなたは車を運転しており、あるタイミングから一時間毎に車の走行速度を計測していたとします。そして、その計測結果が下記であったとします（単位は km / h）。

この時、速度の計測を始めてから２時間後〜８時間後の間で車が走行した距離は何 km であったと考えられるでしょうか？

ご存じの方も多いかもしれませんが、距離は速度の積分で求めることが出来ますので、このような問題も数値積分から求めることが可能です。また、速度と距離の関係だけでなく、電力から電力量を求めるようなときも積分が利用できます。

そして、上記のような速度が計測された場合、特定の時間に車が走行した距離を計算するプログラムのソースコードは下記となります。

#include <stdio.h>

typedef double (*FUNC)(int);

double func(int x) {
    static double v[10] = {52, 60, 58, 20, 35, 0, 43, 29, 51, 27};  
    return v[x];
}

double integral(int a, int b, int n, FUNC f) {
    int h = (b - a) / n;

    double sum = 0;
    for (int i = 0; i < n; i++) {
        sum += 0.5 * (f(a + i * h) + f(a + ((i + 1) * h))) * h;
    }
    return sum;
}

int main(void) {
    double result = integral(2, 8, 8 - 2, func);
    printf("%lf\n", result); 
}

関数に対する数値積分 で示したソースコードから若干変更していますが、積分結果が得られる原理は先程と同じになります。ただ、今回の場合は分割後の各台形の高さは計測タイミングの間隔である１時間で固定となります。

また、積分対象の関数 func は、単純に引数 x を添字とする配列 v の要素 v[x] となります。

上記のソースコードをコンパイルして実行すれば、下記のような結果が得られるはずです。これが、配列 v を積分区間 [2, 8] で積分した結果となります。

181.500000

このように、数値積分であれば、定期的にサンプリングして生成した数列を積分するようなことも可能となります。

ただし、上記の結果にはやはり誤差が含まれます。

今回の例では１時間毎にしか速度を計測していないため、その間の速度の変化が積分結果に反映されないことになります。

例えば、計測したタイミング以外の時間、実は車は停止していました…、という話であれば、実際の走行距離は先ほど求めた数値積分の結果と大きく異なることになるでしょう。

この誤差も、結局は台形への分割数を増やすことで小さくすることができます。今回の場合は、速度を計測する間隔を１時間ではなく１分や１秒にしてやれば、もっと正確に距離が求まることになります。これは直感的にも理解していただけると思います。

このように、数値積分では、誤差を減らすためには求める領域の分割数が非常に重要となります。ただ、分割数を増やすことは計算時間の増加に繋がりますし、今回の場合であれば計測した速度を記憶するためのメモリも増えることになります。積分の精度と計算時間＆メモリ使用量はトレードオフの関係となりますので、このあたりをバランス良く選択することが数値積分を上手く利用するポイントになると思います。

まあ、これは数値積分だけでなく色んなケースで言えることではあるのですが…。

スポンサーリンク

まとめ

このページでは、C言語における台形公式による数値積分のやり方について説明しました！

プログラミングにおいても積分が活躍することが多くあります。そして、この積分は、不定積分などのような難しいことを行わなくても、単純に元々の関数を数値として評価し、その結果から台形公式を用いた計算で求めることが可能です。

この考え方で、速度から距離や、瞬時電力から使用電力量を計算するようなことも可能です。数値積分だけでなく、台形で近似する考え方は色んな時に活用できますので、是非これらの考え方については覚えておきましょう！

オススメの参考書（PR）

C言語学習中だけど分からないことが多くて挫折しそう...という方には、下記の「スッキリわかるC言語入門」がオススメです！

スッキリわかるC言語入門 第2版

posted with ヨメレバ

中山 清喬/株式会社フレアリンク インプレス 2021年12月14日頃

楽天ブックス

Amazon

Kindle

 

まず学習を進める上で、参考書は２冊持っておくことをオススメします。この理由は下記の２つです。

• 参考書によって、解説の仕方は異なる
• 読み手によって、理解しやすい解説の仕方は異なる

ある人の説明聞いても理解できなかったけど、他の人からちょっと違った観点での説明を聞いて「あー、そういうことね！」って簡単に理解できた経験をお持ちの方も多いのではないでしょうか？

それと同じで、１冊の参考書を読んで理解できない事も、他の参考書とは異なる内容の解説を読むことで理解できる可能性があります。

なので、参考書は２冊持っておいた方が学習時に挫折しにくいというのが私の考えです。

特に上記の「スッキリわかるC言語入門」は、他の参考書とは違った切り口での解説が豊富で、他の参考書で理解できなかった内容に対して違った観点での解説を読むことができ、オススメです。題名の通り「なぜそうなるのか？」がスッキリ理解できるような解説内容にもなっており、C言語入門書としてもかなり分かりやすい参考書だと思います。

もちろんネット等でも色んな観点からの解説を読むことが出来ますので、分からない点は別の人・別の参考書の解説を読んで解決していきましょう！もちろん私のサイトも参考にしていただけると嬉しいです！

入門用のオススメ参考書は下記ページでも紹介していますので、こちらも是非参考にしていただければと思います。

https://daeudaeu.com/c_reference_book/