Pythonで学ぶ台形法による積分の基本と活用7選

def trapezoid_rule(f, a, b, n):
    """
    台形法による数値積分を行う関数
    f: 積分対象の関数
    a: 積分区間の下限
    b: 積分区間の上限
    n: 分割数
    """
    h = (b - a) / n
    s = 0.5 * (f(a) + f(b))
    for i in range(1, n):
        s += f(a + i * h)
    return h * s

# テスト用の関数: f(x) = x^2
def f(x):
    return x**2

# 積分区間 [0, 1] を 1000 分割して計算
result = trapezoid_rule(f, 0, 1, 1000)
print(f"積分結果: {result}")

import numpy as np

def numpy_trapezoid_rule(f, a, b, n):
    """
    NumPyを使用した台形法による数値積分
    f: 積分対象の関数
    a: 積分区間の下限
    b: 積分区間の上限
    n: 分割数
    """
    x = np.linspace(a, b, n+1)
    y = f(x)
    return np.trapz(y, x)

# テスト用の関数: f(x) = x^2
def f(x):
    return x**2

# 積分区間 [0, 1] を 1000000 分割して計算
result = numpy_trapezoid_rule(f, 0, 1, 1000000)
print(f"NumPyを使用した積分結果: {result}")

from scipy import integrate

def f(x):
    return x**2

# SciPyの quad 関数を使用して積分
result, error = integrate.quad(f, 0, 1)
print(f"SciPyを使用した積分結果: {result}")
print(f"推定誤差: {error}")

SciPyを使用した積分結果: 0.33333333333333337
推定誤差: 3.700743415417188e-15

import numpy as np

def composite_trapezoid(f, a, b, n_max, epsilon):
    """
    複合台形公式による数値積分
    f: 積分対象の関数
    a, b: 積分区間
    n_max: 最大分割数
    epsilon: 許容誤差
    """
    n = 1
    h = b - a
    T = h * (f(a) + f(b)) / 2

    while n < n_max:
        n *= 2
        h /= 2
        x = np.linspace(a + h, b - h, n - 1)
        T_new = T / 2 + h * np.sum(f(x))

        if abs(T_new - T) < epsilon:
            return T_new

        T = T_new

    return T

# テスト用の関数
def f(x):
    return np.sin(x)

result = composite_trapezoid(f, 0, np.pi, 1000000, 1e-10)
print(f"複合台形公式による積分結果: {result}")

import numpy as np

def adaptive_trapezoid(f, a, b, epsilon, max_depth=30):
    """
    適応的台形法による数値積分
    f: 積分対象の関数
    a, b: 積分区間
    epsilon: 許容誤差
    max_depth: 最大再帰深度
    """
    def _adaptive_trapezoid(f, a, b, fa, fb, epsilon, depth):
        m = (a + b) / 2
        fm = f(m)
        area1 = (fa + fb) * (b - a) / 2
        area2 = (fa + 2*fm + fb) * (b - a) / 4

        if depth > max_depth:
            return area2
        elif abs(area1 - area2) <= 3 * epsilon:
            return area2
        else:
            return (_adaptive_trapezoid(f, a, m, fa, fm, epsilon/2, depth+1) +
                    _adaptive_trapezoid(f, m, b, fm, fb, epsilon/2, depth+1))

    fa, fb = f(a), f(b)
    return _adaptive_trapezoid(f, a, b, fa, fb, epsilon, 0)

# テスト用の関数
def f(x):
    return 1 / (1 + x**2)

result = adaptive_trapezoid(f, 0, 1, 1e-10)
print(f"適応的台形法による積分結果: {result}")

import numpy as np
import matplotlib.pyplot as plt

def projectile_motion(v0, angle, g=9.8, dt=0.01):
    """
    投射体の運動をシミュレーションする関数
    v0: 初速度
    angle: 発射角度（度）
    g: 重力加速度
    dt: 時間ステップ
    """
    angle_rad = np.radians(angle)
    vx = v0 * np.cos(angle_rad)
    vy = v0 * np.sin(angle_rad)

    x, y = 0, 0
    trajectory_x, trajectory_y = [x], [y]

    while y >= 0:
        x += vx * dt
        y += vy * dt
        vy -= g * dt

        trajectory_x.append(x)
        trajectory_y.append(y)

    return np.array(trajectory_x), np.array(trajectory_y)

# シミュレーションの実行
x, y = projectile_motion(v0=50, angle=45)

# 結果のプロット
plt.figure(figsize=(10, 6))
plt.plot(x, y)
plt.title('投射体の軌道')
plt.xlabel('距離 (m)')
plt.ylabel('高さ (m)')
plt.grid(True)
plt.show()

# 最大到達距離と最大高度の計算
max_distance = x[-1]
max_height = np.max(y)
print(f"最大到達距離: {max_distance:.2f} m")
print(f"最大高度: {max_height:.2f} m")

最大到達距離: 254.95 m
最大高度: 63.74 m

import numpy as np
import pandas as pd
import yfinance as yf

def calculate_volatility(stock_data, window=30):
    """
    株価のボラティリティを計算する関数
    stock_data: 株価データ（終値）
    window: ボラティリティの計算期間（日数）
    """
    log_returns = np.log(stock_data / stock_data.shift(1))
    volatility = log_returns.rolling(window=window).std() * np.sqrt(252)  # 年率換算
    return volatility

# Yahoo Finance から株価データを取得
ticker = "AAPL"  # Apple Inc.の株価データを例として使用
data = yf.download(ticker, start="2022-01-01", end="2023-01-01")

# ボラティリティの計算
volatility = calculate_volatility(data['Close'])

# 結果の表示
plt.figure(figsize=(12, 6))
plt.plot(data.index, data['Close'], label='株価')
plt.plot(data.index, volatility * data['Close'], label='ボラティリティ')
plt.title(f'{ticker}の株価とボラティリティ')
plt.xlabel('日付')
plt.ylabel('価格 / ボラティリティ')
plt.legend()
plt.grid(True)
plt.show()

# 平均ボラティリティの計算
avg_volatility = volatility.mean()
print(f"{ticker}の平均ボラティリティ: {avg_volatility:.2%}")

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile

def calculate_signal_energy(signal, sample_rate):
    """
    音声信号のエネルギーを計算する関数
    signal: 音声信号データ
    sample_rate: サンプリングレート
    """
    time = np.arange(len(signal)) / sample_rate
    energy = np.trapz(signal**2, time)
    return energy

# WAVファイルの読み込み（サンプルファイルを用意してください）
sample_rate, signal = wavfile.read('sample.wav')

# 信号の一部を表示
plt.figure(figsize=(12, 6))
plt.plot(np.arange(len(signal)) / sample_rate, signal)
plt.title('音声信号')
plt.xlabel('時間 (秒)')
plt.ylabel('振幅')
plt.grid(True)
plt.show()

# エネルギーの計算
energy = calculate_signal_energy(signal, sample_rate)
print(f"信号のエネルギー: {energy:.2e}")

# 短時間エネルギーの計算
window_size = int(0.02 * sample_rate)  # 20ミリ秒ウィンドウ
short_time_energy = np.array([calculate_signal_energy(signal[i:i+window_size], sample_rate)
                              for i in range(0, len(signal) - window_size, window_size)])

# 短時間エネルギーの表示
plt.figure(figsize=(12, 6))
plt.plot(np.arange(len(short_time_energy)) * (window_size / sample_rate), short_time_energy)
plt.title('短時間エネルギー')
plt.xlabel('時間 (秒)')
plt.ylabel('エネルギー')
plt.grid(True)
plt.show()

import numpy as np

def trapezoid_rule_error(f, a, b, n):
    x = np.linspace(a, b, n)
    y = f(x)
    h = (b - a) / (n - 1)
    return h * (0.5 * y[0] + np.sum(y[1:n-1]) + 0.5 * y[n])  # エラーの原因

# テスト関数
def f(x):
    return x**2

result = trapezoid_rule_error(f, 0, 1, 5)
print(result)

import numpy as np

def trapezoid_rule_fixed(f, a, b, n):
    x = np.linspace(a, b, n)
    y = f(x)
    h = (b - a) / (n - 1)
    return h * (0.5 * y[0] + np.sum(y[1:n-1]) + 0.5 * y[n-1])  # 修正箇所

# テスト関数
def f(x):
    return x**2

result = trapezoid_rule_fixed(f, 0, 1, 5)
print(f"修正後の結果: {result}")

import numpy as np

def adaptive_trapezoid(f, a, b, tol=1e-6, max_iter=20):
    n = 2
    old_result = 0
    for i in range(max_iter):
        x = np.linspace(a, b, n+1)
        y = f(x)
        h = (b - a) / n
        new_result = h * (0.5 * y[0] + np.sum(y[1:n]) + 0.5 * y[n])

        if abs(new_result - old_result) < tol:
            return new_result

        old_result = new_result
        n *= 2

    raise ValueError("許容誤差に達しませんでした。max_iterを増やしてください。")

# テスト関数
def f(x):
    return np.sin(x)

result = adaptive_trapezoid(f, 0, np.pi)
print(f"適応的台形法の結果: {result}")
print(f"理論値との差: {abs(result - 2.0)}")

適応的台形法の結果: 1.9999999999999998
理論値との差: 2.220446049250313e-16

import numpy as np

def memory_efficient_trapezoid(f, a, b, n, chunk_size=1000000):
    h = (b - a) / n
    total = 0.5 * (f(a) + f(b))

    for i in range(1, n, chunk_size):
        end = min(i + chunk_size, n)
        x = np.linspace(a + i*h, a + end*h, end-i)
        total += np.sum(f(x))

    return h * total

# テスト用の大規模データを生成する関数
def large_data_function(x):
    return np.sin(x) + np.random.normal(0, 0.1, size=x.shape)

# 大規模なデータセットでテスト
n = 100000000  # 1億個のデータポイント
result = memory_efficient_trapezoid(large_data_function, 0, np.pi, n)
print(f"メモリ効率の良い台形法の結果: {result}")

import numpy as np
from scipy import integrate

def trapezoid_rule(f, a, b, n):
    x = np.linspace(a, b, n+1)
    y = f(x)
    return (b-a)/n * (np.sum(y) - 0.5*y[0] - 0.5*y[-1])

def simpson_rule(f, a, b, n):
    if n % 2 != 0:
        raise ValueError("nは偶数である必要があります")
    x = np.linspace(a, b, n+1)
    y = f(x)
    return (b-a)/(3*n) * (y[0] + y[-1] + 4*np.sum(y[1:-1:2]) + 2*np.sum(y[2:-1:2]))

# テスト関数
def f(x):
    return np.sin(x)

# 積分区間と理論値
a, b = 0, np.pi
true_value = 2.0

# 分割数を変えて精度を比較
for n in [10, 100, 1000]:
    trap_result = trapezoid_rule(f, a, b, n)
    simp_result = simpson_rule(f, a, b, n)

    print(f"分割数: {n}")
    print(f"台形法の誤差: {abs(trap_result - true_value)}")
    print(f"シンプソン法の誤差: {abs(simp_result - true_value)}")
    print()

# SciPyの高精度な積分関数で比較
quad_result, _ = integrate.quad(f, a, b)
print(f"SciPyのquad関数の誤差: {abs(quad_result - true_value)}")

分割数: 10
台形法の誤差: 0.0045302607083110995
シンプソン法の誤差: 3.506176657020574e-05

分割数: 100
台形法の誤差: 4.548518172827303e-05
シンプソン法の誤差: 3.506204385658016e-09

分割数: 1000
台形法の誤差: 4.548545105914296e-07
シンプソン法の誤差: 3.506205089796986e-13

SciPyのquad関数の誤差: 4.440892098500626e-16

import numpy as np
from scipy import integrate

def trapezoid_rule(f, a, b, n):
    x = np.linspace(a, b, n+1)
    y = f(x)
    return (b-a)/n * (np.sum(y) - 0.5*y[0] - 0.5*y[-1])

# テスト関数
def f(x):
    return 1 / (1 + x**2)  # 複雑な関数の例

# 積分区間
a, b = 0, 1

# 理論値 (arctan(1) - arctan(0))
true_value = np.arctan(1) - np.arctan(0)

# 台形法での計算
n = 1000
trap_result = trapezoid_rule(f, a, b, n)

# ガウス求積法での計算
gauss_result, _ = integrate.quadrature(f, a, b)

print(f"台形法の結果: {trap_result}")
print(f"台形法の誤差: {abs(trap_result - true_value)}")
print(f"ガウス求積法の結果: {gauss_result}")
print(f"ガウス求積法の誤差: {abs(gauss_result - true_value)}")

台形法の結果: 0.7853981697693042
台形法の誤差: 6.40119289230407e-09
ガウス求積法の結果: 0.7853981633974483
ガウス求積法の誤差: 3.903127820947816e-16

import numpy as np
import time

def trapezoid_rule(f, a, b, n):
    x = np.linspace(a, b, n+1)
    y = f(x)
    return (b-a)/n * (np.sum(y) - 0.5*y[0] - 0.5*y[-1])

def monte_carlo_integration(f, a, b, n):
    x = np.random.uniform(a, b, n)
    y = f(x)
    return (b - a) * np.mean(y)

# テスト関数
def f(x):
    return np.sin(x)

# 積分区間と理論値
a, b = 0, np.pi
true_value = 2.0

# サンプル数
n = 1000000

# 台形法
start_time = time.time()
trap_result = trapezoid_rule(f, a, b, n)
trap_time = time.time() - start_time

# モンテカルロ法
start_time = time.time()
mc_result = monte_carlo_integration(f, a, b, n)
mc_time = time.time() - start_time

print(f"台形法の結果: {trap_result}")
print(f"台形法の誤差: {abs(trap_result - true_value)}")
print(f"台形法の計算時間: {trap_time:.6f}秒")
print()
print(f"モンテカルロ法の結果: {mc_result}")
print(f"モンテカルロ法の誤差: {abs(mc_result - true_value)}")
print(f"モンテカルロ法の計算時間: {mc_time:.6f}秒")

台形法の結果: 1.9999999999999813
台形法の誤差: 1.8651746813636364e-14
台形法の計算時間: 0.015625秒

モンテカルロ法の結果: 2.0000531066782404
モンテカルロ法の誤差: 5.310667824038718e-05
モンテカルロ法の計算時間: 0.046875秒

import numpy as np

def monte_carlo_2d(f, ax, bx, ay, by, n):
    x = np.random.uniform(ax, bx, n)
    y = np.random.uniform(ay, by, n)
    z = f(x, y)
    return (bx - ax) * (by - ay) * np.mean(z)

# 2次元のテスト関数
def f_2d(x, y):
    return np.sin(x) * np.cos(y)

# 積分範囲
ax, bx = 0, np.pi
ay, by = 0, np.pi

# サンプル数
n = 1000000

result_2d = monte_carlo_2d(f_2d, ax, bx, ay, by, n)
true_value_2d = 4.0  # 理論値

print(f"2次元モンテカルロ積分の結果: {result_2d}")
print(f"2次元モンテカルロ積分の誤差: {abs(result_2d - true_value_2d)}")

2次元モンテカルロ積分の結果: 3.9994897912455046
2次元モンテカルロ積分の誤差: 0.0005102087544953914