Juliaで数値積分: 台形公式

Juliaプログラミングに入る前に、まずは数値積分の基本的な概念を探りましょう。数値積分は、数値計算を使用して関数の積分を評価する手法を含みます。次の一次元積分を考えてみましょう： \[ S=\int_{a}^{b} dx \, f(x). \tag{1} \] ここで、積分区間の \(N\)-分割を想像してみましょう（\(a \le x \le b\)）。各分割の幅 \(\Delta x\) は \(\Delta x = (b - a)/N\) と定義されます。積分 \(S\) は各分割区間ごとの積分の合計として表現できます： \begin{align} S &= \int_{a}^{b} dx \, f(x) \\ &= \int_{a+\Delta x}^{a} dx \, f(x) + \int_{a+2\Delta x}^{a+\Delta x} dx \, f(x) + \int_{a+3\Delta x}^{a+2\Delta x} dx \, f(x) + \cdots + \int_{b}^{a+(N-1)\Delta x} dx \, f(x) \\ &= \sum_{j=0}^{N-1} \int_{a + j\Delta x}^{a + (j+1)\Delta x} dx \, f(x) \tag{2} \end{align}

分割幅 \(\Delta x\) が十分に小さいと仮定すると、被積分関数 \(f(x)\) は各小さな積分区間の端の値の平均として近似できます。言い換えれば、式 (2) は次のように近似できます： \begin{align} S &= \sum_{j=0}^{N-1} \int_{a + j\Delta x}^{a + (j+1)\Delta x} f(x) dx \\ & \approx \sum_{j=0}^{N-1} \int_{a + j\Delta x}^{a + (j+1)\Delta x} \frac{f \left ( a + (j+1)\Delta x \right )+f\left (a + j\Delta x \right )}{2} dx \\ &= \sum_{j=0}^{N-1}\frac{f \left ( a + (j+1)\Delta x \right )+f\left (a + j\Delta x \right )}{2} \int_{a + j\Delta x}^{a + (j+1)\Delta x} \, dx \\ &= \sum_{j=0}^{N-1} \frac{f \left ( a + (j+1)\Delta x \right )+f\left (a + j\Delta x \right )}{2} \, \Delta x \\ &= \frac{f(a)}{2}\Delta x + \left [ \sum_{j=1}^{N-1} f(a + j\Delta x) \, \Delta x \right ] + \frac{f(b)}{2}\Delta x \tag{3} \end{align}

式（3）で示される式は一般的に台形則として知られており、これは台形の面積を合計することで積分値を近似するためのものです。区間の分割数が大きい場合（\(N\rightarrow \infty\)、\(\Delta x \rightarrow 0\)）、式（3）の台形則は正確な積分値 \(S\) に収束します。

練習として、次の積分を台形則を用いて計算してみましょう： \begin{align} S = \int_{0}^{1} dx \, \frac{4}{1+x^2} = \pi . \end{align} 以下に示すのは、この積分を数値的に評価するJuliaコードです。数値積分とJuliaプログラミングに慣れるために、提供されたサンプルコードを元に自分自身のコードを作成してみてください。

(lesson2_1.jl)

# define your function
function my_func(x)
    return 4.0 / (1.0 + x^2)
end

# integrate the function from a to b with n points
function integrate(f, a, b, n)
    h = (b - a) / n

    s = (f(a) + f(b)) / 2.0
    for i in 1:n-1
        x = a + i * h
        s += f(x)
    end

    s *= h
    return s
end

a = 0.0
b = 1.0
n = 64

num_integral = integrate(my_func, a, b, n)

println("num. integral = ", num_integral)
println("pi            = ", 1.0*π)

提供されたJuliaコードは、2つの関数を定義し、特定の数学的関数の数値積分の近似を実行します。

最初の関数 my_func(x) は、数値入力 x を受け取り、\(\frac{4.0}{1.0+x^2}\) の値を計算します。
2番目の関数 integrate(f, a, b, n) は、指定された関数 f の区間 \([a, b]\) での数値積分を、点数 n を用いた台形則を使用して計算します。各台形の幅を変数 h として初期化し、区間端点での関数の平均高さを計算し、内部の点を反復処理して、それらの寄与を合計面積に加えます。最終的に、各台形の幅で結果を乗じて近似積分値を得ます。

これらの関数を定義した後、コードは区間 \([a, b]\)（0.0から1.0）および点数 n（64）の具体的な値を設定します。その後、関数 my_func と指定された区間および点数で integrate 関数を呼び出し、その結果を変数 num_integral に格納します。

最後に、コードは数値積分結果（num_integral）を表示し、それをπの値と比較して両方の値を表示します。

Lesson 2: 数値積分