【C言語】積分するプログラムの書き方

2023-09-132024-07-22

数値積分とは、関数の面積を数値的に求める手法です。

具体的には、台形公式とシンプソンの公式という2つの代表的な手法を使って積分を行います。

さらに、これらの手法を用いたC言語のプログラムの実装方法や、実際のデータを使った応用例についても解説します。

目次から探す

数値積分の手法

数値積分は、関数の積分を数値的に求める手法です。

数値積分の手法にはいくつかの種類がありますが、ここでは代表的な「台形公式」と「シンプソンの公式」について解説します。

台形公式

台形公式の原理

台形公式は、積分区間を小さな台形に分割し、それぞれの台形の面積を合計することで積分値を求める手法です。

具体的には、関数 ( f(x) ) を積分区間 ([a, b]) で ( n ) 等分し、各区間の端点での関数値を用いて台形の面積を計算します。

台形公式の数式は以下のようになります：

ここで、です。

台形公式の利点と欠点

利点:

実装が簡単で理解しやすい。
任意の関数に対して適用可能。

欠点:

精度が低い場合がある。
特に関数が急激に変化する場合には誤差が大きくなる。
分割数 ( n ) を増やすと計算量が増加する。

シンプソンの公式

シンプソンの公式の原理

シンプソンの公式は、積分区間を小さなパラボラに分割し、それぞれのパラボラの面積を合計することで積分値を求める手法です。

具体的には、関数 ( f(x) ) を積分区間 ([a, b]) で ( n ) 等分し、各区間の端点と中点での関数値を用いてパラボラの面積を計算します。

シンプソンの公式の数式は以下のようになります：

ここで、です。

シンプソンの公式の利点と欠点

利点:

台形公式よりも高い精度を持つ。
分割数 ( n ) が少なくても比較的高い精度を保つ。

欠点:

実装が台形公式よりも複雑。
分割数 ( n ) は偶数でなければならない。

以上が、数値積分の代表的な手法である「台形公式」と「シンプソンの公式」です。

次に、これらの手法を用いたC言語のプログラムの実装方法について詳しく見ていきます。

台形公式を用いた積分プログラム

台形公式の実装

台形公式は、関数のグラフを台形で近似することで積分を行う手法です。

具体的には、積分区間を小さな区間に分割し、それぞれの区間で台形の面積を計算して合計します。

関数の定義

まず、積分する関数と台形公式を用いた積分を行う関数を定義します。

ここでは、例として関数 f(x) = x^2 を積分します。

#include <stdio.h>
// 積分する関数 f(x) = x^2
double f(double x) {
    return x * x;
}
// 台形公式を用いた積分関数
double trapezoidal(double (*func)(double), double a, double b, int n) {
    double h = (b - a) / n; // 区間の幅
    double sum = 0.5 * (func(a) + func(b)); // 端点の処理
    for (int i = 1; i < n; i++) {
        sum += func(a + i * h);
    }
    return sum * h;
}

メイン関数での呼び出し

次に、メイン関数で上記の関数を呼び出して積分を実行します。

積分区間 [a, b] と分割数 n を指定します。

int main() {
    double a = 0.0; // 積分区間の開始点
    double b = 1.0; // 積分区間の終了点
    int n = 100; // 分割数
    double result = trapezoidal(f, a, b, n);
    printf("積分結果: %f\n", result);
    return 0;
}

サンプルコード

上記の内容をまとめた完全なサンプルコードは以下の通りです。

#include <stdio.h>
// 積分する関数 f(x) = x^2
double f(double x) {
    return x * x;
}
// 台形公式を用いた積分関数
double trapezoidal(double (*func)(double), double a, double b, int n) {
    double h = (b - a) / n; // 区間の幅
    double sum = 0.5 * (func(a) + func(b)); // 端点の処理
    for (int i = 1; i < n; i++) {
        sum += func(a + i * h);
    }
    return sum * h;
}
int main() {
    double a = 0.0; // 積分区間の開始点
    double b = 1.0; // 積分区間の終了点
    int n = 100; // 分割数
    double result = trapezoidal(f, a, b, n);
    printf("積分結果: %f\n", result);
    return 0;
}

実行結果の確認

上記のプログラムを実行すると、以下のような結果が得られます。

積分結果: 0.333350

この結果は、関数 f(x) = x^2 を区間 [0, 1] で積分した値に近いことがわかります。

理論的な正確な値は 1/3 (約 0.333333) です。

分割数 n を増やすことで、より精度の高い結果が得られます。

シンプソンの公式を用いた積分プログラム

シンプソンの公式は、数値積分の手法の一つで、台形公式よりも精度が高いとされています。

ここでは、シンプソンの公式を用いた積分プログラムの実装方法について解説します。

シンプソンの公式の実装

シンプソンの公式を用いた積分プログラムを実装するためには、まず関数の定義を行い、その後メイン関数で呼び出す形にします。

関数の定義

シンプソンの公式を用いた積分を行う関数を定義します。

この関数は、積分区間の開始点、終了点、および分割数を引数として受け取ります。

#include <stdio.h>
// 積分する関数の定義
double f(double x) {
    return x * x; // 例として f(x) = x^2 を使用
}
// シンプソンの公式を用いた積分関数の定義
double simpsons_rule(double a, double b, int n) {
    double h = (b - a) / n;
    double sum = f(a) + f(b);
    for (int i = 1; i < n; i += 2) {
        sum += 4 * f(a + i * h);
    }
    for (int i = 2; i < n - 1; i += 2) {
        sum += 2 * f(a + i * h);
    }
    return (h / 3) * sum;
}

メイン関数での呼び出し

次に、メイン関数でシンプソンの公式を用いた積分関数を呼び出します。

ここでは、積分区間を0から1、分割数を100とします。

int main() {
    double a = 0.0; // 積分区間の開始点
    double b = 1.0; // 積分区間の終了点
    int n = 100;    // 分割数
    double result = simpsons_rule(a, b, n);
    printf("積分結果: %f\n", result);
    return 0;
}

サンプルコード

上記の関数定義とメイン関数を組み合わせた完全なサンプルコードは以下の通りです。

#include <stdio.h>
// 積分する関数の定義
double f(double x) {
    return x * x; // 例として f(x) = x^2 を使用
}
// シンプソンの公式を用いた積分関数の定義
double simpsons_rule(double a, double b, int n) {
    double h = (b - a) / n;
    double sum = f(a) + f(b);
    for (int i = 1; i < n; i += 2) {
        sum += 4 * f(a + i * h);
    }
    for (int i = 2; i < n - 1; i += 2) {
        sum += 2 * f(a + i * h);
    }
    return (h / 3) * sum;
}
int main() {
    double a = 0.0; // 積分区間の開始点
    double b = 1.0; // 積分区間の終了点
    int n = 100;    // 分割数
    double result = simpsons_rule(a, b, n);
    printf("積分結果: %f\n", result);
    return 0;
}

実行結果の確認

上記のサンプルコードを実行すると、以下のような結果が得られます。

積分結果: 0.333333

この結果は、関数 f(x) = x^2 を区間 [0, 1] で積分した結果が 1/3 であることを示しています。

シンプソンの公式を用いることで、非常に高い精度で数値積分を行うことができることがわかります。

以上で、シンプソンの公式を用いた積分プログラムの実装方法についての解説を終わります。

次に、応用例について見ていきましょう。

応用例

数値積分の基本的な手法を理解したところで、次に応用例について見ていきましょう。

ここでは、複数の関数の積分、実際のデータを用いた積分、そしてグラフ表示による結果の可視化について解説します。

複数関数の積分

複数の関数を積分する場合、基本的な手法は単一の関数を積分する場合と同じです。

ただし、異なる関数をそれぞれ積分するために、関数の定義や呼び出し方法を工夫する必要があります。

サンプルコード

以下に、2つの関数 ( f(x) = x^2 ) と ( g(x) = \sin(x) ) を台形公式で積分するサンプルコードを示します。

#include <stdio.h>
#include <math.h>
// 関数 f(x) = x^2
double f(double x) {
    return x * x;
}
// 関数 g(x) = sin(x)
double g(double x) {
    return sin(x);
}
// 台形公式による積分
double trapezoidal(double (*func)(double), double a, double b, int n) {
    double h = (b - a) / n;
    double sum = (func(a) + func(b)) / 2.0;
    for (int i = 1; i < n; i++) {
        sum += func(a + i * h);
    }
    return sum * h;
}
int main() {
    double a = 0.0, b = 1.0;
    int n = 1000;
    double integral_f = trapezoidal(f, a, b, n);
    double integral_g = trapezoidal(g, a, b, n);
    printf("Integral of f(x) = x^2 from %f to %f is %f\n", a, b, integral_f);
    printf("Integral of g(x) = sin(x) from %f to %f is %f\n", a, b, integral_g);
    return 0;
}

実行結果の確認

このプログラムを実行すると、以下のような結果が得られます。

Integral of f(x) = x^2 from 0.000000 to 1.000000 is 0.333333
Integral of g(x) = sin(x) from 0.000000 to 1.000000 is 0.459698

実際のデータを用いた積分

実際のデータを用いた積分では、データポイントを数値積分の入力として使用します。

例えば、実験データやセンサーデータを積分することで、累積量を求めることができます。

サンプルコード

以下に、実際のデータを用いた積分のサンプルコードを示します。

ここでは、データポイントが配列として与えられていると仮定します。

#include <stdio.h>
// 実際のデータポイント
double data[] = {0.0, 0.1, 0.4, 0.9, 1.6, 2.5, 3.6, 4.9, 6.4, 8.1};
int data_size = sizeof(data) / sizeof(data[0]);
// 台形公式による積分
double trapezoidal_data(double *data, int size, double h) {
    double sum = (data[0] + data[size - 1]) / 2.0;
    for (int i = 1; i < size - 1; i++) {
        sum += data[i];
    }
    return sum * h;
}
int main() {
    double h = 1.0; // データポイント間の間隔
    double integral = trapezoidal_data(data, data_size, h);
    printf("Integral of the data points is %f\n", integral);
    return 0;
}

実行結果の確認

このプログラムを実行すると、以下のような結果が得られます。

Integral of the data points is 24.500000

グラフ表示による結果の可視化

積分結果をグラフで表示することで、視覚的に理解しやすくなります。

C言語自体にはグラフ表示の機能はありませんが、外部ライブラリやツールを使用することで実現可能です。

C言語でグラフを出力するのは少々難しい(別のライブラリの知識が必要である)ためここでは、Pythonのmatplotlibを使用してグラフを表示する方法を紹介します。

サンプルコード

まず、C言語で積分結果をファイルに出力します。

#include <stdio.h>
#include <math.h>
// 関数 f(x) = x^2
double f(double x) {
    return x * x;
}
// 台形公式による積分
double trapezoidal(double (*func)(double), double a, double b, int n) {
    double h = (b - a) / n;
    double sum = (func(a) + func(b)) / 2.0;
    for (int i = 1; i < n; i++) {
        sum += func(a + i * h);
    }
    return sum * h;
}
int main() {
    double a = 0.0, b = 1.0;
    int n = 1000;
    double integral = trapezoidal(f, a, b, n);
    // 結果をファイルに出力
    FILE *file = fopen("integral_result.txt", "w");
    fprintf(file, "Integral of f(x) = x^2 from %f to %f is %f\n", a, b, integral);
    fclose(file);
    return 0;
}

次に、Pythonでファイルを読み込み、グラフを表示します。

import matplotlib.pyplot as plt
import numpy as np
# C言語で出力した結果を読み込む
with open("integral_result.txt", "r") as file:
    result = file.readline().strip()
# グラフを表示
x = np.linspace(0, 1, 100)
y = x**2
plt.plot(x, y, label='f(x) = x^2')
plt.fill_between(x, y, alpha=0.2)
plt.title(result)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.legend()
plt.show()