科学技術計算
はじめに
[編集]本書は、科学技術計算に携わる研究者、技術者、学生のための実践的なハンドブックです。基礎的な数学から最新の機械学習手法まで、幅広い内容を網羅しています。各章には具体的な計算例とPythonによるプログラミング例を収録し、実務での即戦力となる知識の習得を目指しています。
第1章 基礎数学
[編集]1.1 代数学の基礎
[編集]数学的な記述の基礎となる代数学について解説します。特に、科学技術計算で頻繁に現れる多項式、行列、ベクトルの取り扱いに重点を置きます。
1.1.1 集合と写像
[編集]集合論の基本概念から始めましょう。集合とは、ものの集まりを表す最も基本的な概念です。
- 定義
- 集合
集合とは、明確な基準によって、ものの属する・属さないが判定できる、ものの集まりです。
例えば、以下のような集合が考えられます：
- 自然数の集合 N = {1, 2, 3, ...}
- 整数の集合 Z = {..., -2, -1, 0, 1, 2, ...}
- 実数の集合 R
- 集合の演算
集合A, Bに対して、以下の基本的な演算が定義されます：
- 和集合 A ∪ B：AまたはBのいずれかに属する要素の集合
- 積集合 A ∩ B：AとBの両方に属する要素の集合
- 差集合 A \ B：Aに属し、Bに属さない要素の集合
1.1.2 群・環・体
[編集]代数的構造の基本概念である群、環、体について説明します。
- 群の定義
集合Gと二項演算・が与えられているとき、以下の4つの条件を満たすとき、(G, ・)を群といいます：
- 結合法則：任意の a, b, c ∈ G に対して (a・b)・c = a・(b・c)
- 単位元の存在：e ∈ G が存在して、任意の a ∈ G に対して a・e = e・a = a
- 逆元の存在：任意の a ∈ G に対して、a・b = b・a = e となる b ∈ G が存在
- 閉じている：任意の a, b ∈ G に対して a・b ∈ G
実践例：
import numpy as np # 回転行列の群の例 def rotation_matrix(theta): return np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]) # 2つの回転の合成 theta1 = np.pi/4 # 45度 theta2 = np.pi/3 # 60度 R1 = rotation_matrix(theta1) R2 = rotation_matrix(theta2) # 結合法則の確認 R3 = rotation_matrix(theta1 + theta2) print(np.allclose(R1 @ R2, R3)) # True
1.2 微分積分学
[編集]1.2.1 微分の基礎
[編集]微分は、関数の局所的な変化率を表す重要な概念です。科学技術計算において、物理現象の記述や最適化問題の解法に不可欠です。
- 導関数の定義
関数 の における導関数 は以下のように定義されます：
- 偏微分
多変数関数 f(x₁, x₂, ..., xₙ) に対する xᵢ についての偏微分は、他の変数を定数とみなして行う微分です：
実践例：
import numpy as np from scipy.misc import derivative # 数値微分の例 def f(x): return x**2 * np.sin(x) # scipy.misc.derivativeを使用した数値微分 x = 2.0 dx = derivative(f, x, dx=1e-6) print(f"f'({x}) ≈ {dx:.6f}") # 中心差分法による独自実装 def numerical_derivative(f, x, h=1e-6): return (f(x + h) - f(x - h)) / (2 * h) dx_custom = numerical_derivative(f, x) print(f"独自実装による f'({x}) ≈ {dx_custom:.6f}")
1.2.2 積分の基礎
[編集]積分は、関数の累積的な効果や面積、体積を計算する際に用いられる基本的な演算です。
- 定積分の定義
関数 f(x) の区間 [a,b] における定積分は以下のように定義されます：
ここで、であり、xᵢ は各小区間の代表点です。
- 数値積分の実装
def trapezoidal_integration(f, a, b, n=1000): """台形則による数値積分""" x = np.linspace(a, b, n) dx = (b - a) / (n - 1) return dx * (np.sum(f(x)) - (f(a) + f(b))/2) # 例：sin(x)の0からπまでの積分（結果は2になるはず） def g(x): return np.sin(x) result = trapezoidal_integration(g, 0, np.pi) print(f"∫[0→π] sin(x)dx ≈ {result:.6f}") # より高度な方法：scipy.integrateの使用 from scipy import integrate result_scipy, error = integrate.quad(g, 0, np.pi) print(f"scipy使用: ∫[0→π] sin(x)dx ≈ {result_scipy:.6f}")
1.2.3 応用例：最適化問題
[編集]微分積分学の重要な応用例として、最適化問題を考えましょう。特に、勾配降下法は機械学習でも頻繁に使用される手法です。
def gradient_descent(f, df, x0, learning_rate=0.1, max_iter=100, tol=1e-6): """ 勾配降下法による最適化 f: 目的関数 df: fの導関数 x0: 初期値 """ x = x0 history = [x] for i in range(max_iter): gradient = df(x) x_new = x - learning_rate * gradient if abs(f(x_new) - f(x)) < tol: break x = x_new history.append(x) return x, history # 例：x²の最小値を求める def f(x): return x**2 def df(x): return 2*x x_min, history = gradient_descent(f, df, x0=2.0) print(f"x²の最小値は x = {x_min:.6f} で達成されます")
1.3 線形代数
[編集]1.3.1 ベクトルと行列の基礎
[編集]線形代数は、科学技術計算の中で最も頻繁に使用される数学的道具の一つです。特に、大規模なデータ処理や機械学習において重要な役割を果たします。
1.3.1 ベクトルと行列の基本演算
[編集]ベクトルと行列は、多次元データを扱う上で基本となる数学的オブジェクトです。
- ベクトルの基本演算
- 内積：a·b = Σᵢ aᵢbᵢ
- ノルム：||a|| = √(Σᵢ aᵢ²)
- 外積（3次元の場合）：a×b = (a₂b₃-a₃b₂, a₃b₁-a₁b₃, a₁b₂-a₂b₁)
import numpy as np # ベクトル演算の例 a = np.array([1, 2, 3]) b = np.array([4, 5, 6]) # 内積 dot_product = np.dot(a, b) print(f"内積: {dot_product}") # ノルム（L2ノルム） norm_a = np.linalg.norm(a) print(f"ベクトルaのノルム: {norm_a:.4f}") # 外積 cross_product = np.cross(a, b) print(f"外積: {cross_product}")
- 行列の基本演算
- 行列の和差：(A±B)ᵢⱼ = Aᵢⱼ±Bᵢⱼ
- 行列の積：(AB)ᵢⱼ = Σₖ AᵢₖBₖⱼ
- 転置：(Aᵀ)ᵢⱼ = Aⱼᵢ
# 行列演算の例 A = np.array([[1, 2], [3, 4]]) B = np.array([[5, 6], [7, 8]]) # 行列の和 C = A + B print("行列の和:\n", C) # 行列の積 D = np.matmul(A, B) # または A @ B print("行列の積:\n", D) # 転置 AT = A.T print("Aの転置:\n", AT)
1.3.2 線形方程式系の解法
[編集]線形方程式系 の解法は、科学技術計算における基本的な問題の一つです。
- ガウス消去法
def gaussian_elimination(A, b): """ ガウス消去法による線形方程式系の解法 A: n×n係数行列 b: n次元ベクトル """ n = len(A) # 拡大係数行列を作成 Ab = np.concatenate((A, b.reshape(n, 1)), axis=1) # 前進消去 for i in range(n): # ピボット選択 pivot = abs(Ab[i:, i]).argmax() + i if pivot != i: Ab[i], Ab[pivot] = Ab[pivot].copy(), Ab[i].copy() for j in range(i + 1, n): factor = Ab[j, i] / Ab[i, i] Ab[j] -= factor * Ab[i] # 後退代入 x = np.zeros(n) for i in range(n-1, -1, -1): x[i] = (Ab[i, -1] - np.dot(Ab[i, i+1:n], x[i+1:])) / Ab[i, i] return x # 例題 A = np.array([[2, 1, -1], [-3, -1, 2], [-2, 1, 2]]) b = np.array([8, -11, -3]) x = gaussian_elimination(A, b) print("解:", x) # 検証：numpy.linalgを使用 x_numpy = np.linalg.solve(A, b) print("numpy.linalg.solveによる解:", x_numpy)
1.3.3 固有値と固有ベクトル
[編集]固有値問題は、振動解析や主成分分析など、多くの応用場面で現れます。
def power_method(A, max_iter=1000, tol=1e-10): """ べき乗法による最大固有値と対応する固有ベクトルの計算 """ n = len(A) x = np.random.rand(n) x = x / np.linalg.norm(x) for _ in range(max_iter): x_new = A @ x x_new = x_new / np.linalg.norm(x_new) if np.allclose(x, x_new, rtol=tol): break x = x_new eigenvalue = (x.T @ A @ x) / (x.T @ x) return eigenvalue, x # 例：対称行列の固有値・固有ベクトル A = np.array([[4, -2], [-2, 5]]) # べき乗法による計算 eval_power, evec_power = power_method(A) print("べき乗法による最大固有値:", eval_power) print("対応する固有ベクトル:", evec_power) # numpy.linalgによる計算（検証用） evals, evecs = np.linalg.eig(A) print("\nnumpy.linalg.eigによる固有値:", evals) print("対応する固有ベクトル:", evecs)
1.3.4 応用例：画像の主成分分析
[編集]線形代数の重要な応用例として、画像の次元削減を行う主成分分析を実装してみましょう。
def pca_compression(image, n_components): """ PCAによる画像圧縮 image: グレースケール画像（2次元配列） n_components: 使用する主成分の数 """ # 平均を引く mean = np.mean(image, axis=1, keepdims=True) centered = image - mean # 共分散行列の計算 cov = centered @ centered.T / centered.shape[1] # 固有値分解 eigenvalues, eigenvectors = np.linalg.eigh(cov) idx = np.argsort(eigenvalues)[::-1] eigenvalues = eigenvalues[idx] eigenvectors = eigenvectors[:, idx] # 次元削減 components = eigenvectors[:, :n_components] projected = components.T @ centered reconstructed = components @ projected + mean return reconstructed # 例：ランダムなパターンの圧縮 np.random.seed(42) original = np.random.rand(100, 1000) # 100×1000のランダム画像 compressed = pca_compression(original, n_components=10) # 圧縮率と再構成誤差の計算 compression_ratio = original.size / (compressed.size + compressed.shape[1]) error = np.mean((original - compressed) ** 2) print(f"圧縮率: {compression_ratio:.2f}") print(f"平均二乗誤差: {error:.6f}")
1.4 複素解析
[編集]複素解析は、科学技術計算において信号処理、制御理論、量子力学など、様々な分野で重要な役割を果たします。
1.4.1 複素数の基礎
[編集]複素数 z = x + yi は実部 x と虚部 y からなり、i は虚数単位（i² = -1）を表します。
- 複素数の基本演算
import numpy as np import cmath # 複素数の生成と基本演算 z1 = 2 + 3j z2 = 1 - 2j # 四則演算 print(f"和: {z1 + z2}") print(f"差: {z1 - z2}") print(f"積: {z1 * z2}") print(f"商: {z1 / z2}") # 極形式への変換 r = abs(z1) # 絶対値（モジュラス） theta = cmath.phase(z1) # 偏角（アーギュメント） print(f"絶対値: {r}") print(f"偏角: {theta} rad") # 極形式から直交形式への変換 z3 = r * (np.cos(theta) + 1j * np.sin(theta)) print(f"極形式から再構成: {z3}")
1.4.2 複素関数と解析性
[編集]複素関数 f(z) が解析的であるとは、その導関数が存在し、連続であることを意味します。
- コーシー・リーマンの方程式
f(z) = u(x,y) + iv(x,y) が解析的であるための必要条件：
def check_analytic(u, v, x, y, h=1e-6): """ コーシー・リーマンの方程式を数値的に確認する """ # 偏微分を数値的に計算 du_dx = (u(x + h, y) - u(x - h, y)) / (2 * h) du_dy = (u(x, y + h) - u(x, y - h)) / (2 * h) dv_dx = (v(x + h, y) - v(x - h, y)) / (2 * h) dv_dy = (v(x, y + h) - v(x, y - h)) / (2 * h) print(f"∂u/∂x ≈ {du_dx:.6f}, ∂v/∂y ≈ {dv_dy:.6f}") print(f"∂u/∂y ≈ {du_dy:.6f}, -∂v/∂x ≈ {-dv_dx:.6f}") return np.allclose(du_dx, dv_dy) and np.allclose(du_dy, -dv_dx) # 例：f(z) = z² = (x + yi)² の場合 def u(x, y): return x'''2 - y'''2 # 実部 def v(x, y): return 2*x*y # 虚部 # z = 1 + i での解析性を確認 is_analytic = check_analytic(u, v, 1.0, 1.0) print(f"関数は解析的か: {is_analytic}")
1.4.3 複素積分
[編集]複素積分は閉曲線に沿った積分として定義され、留数定理は複素積分を計算する強力な道具となります。
def residue_theorem_example(): """ 留数定理を用いた積分の例： ∮ (1/(z²+1)) dz を単位円周上で計算 """ # 極点: z = ±i # 単位円内の極: z = i # 留数: 1/(2i) import scipy.integrate as integrate def integrand(theta): # z = e^(iθ) とパラメータ化 z = np.exp(1j * theta) return 1 / (z**2 + 1) * 1j * z # 数値積分（0 から 2π） result_numerical, _ = integrate.quad( lambda t: np.real(integrand(t)), 0, 2*np.pi ) # 留数定理による結果 result_residue = 2 * np.pi * 1j * (1/(2*1j)) print(f"数値積分結果: {result_numerical:.6f}") print(f"留数定理による結果: {result_residue.real:.6f}") residue_theorem_example()
1.4.4 応用例：FFTを用いた信号処理
[編集]複素解析の重要な応用例として、高速フーリエ変換（FFT）を用いた信号処理を実装してみましょう。
def signal_processing_example(): """ FFTを用いた信号のフィルタリング例 """ # サンプリングパラメータ t = np.linspace(0, 1, 1000) dt = t[1] - t[0] freq = np.fft.fftfreq(len(t), dt) # 元の信号（2つの周波数成分） signal = np.sin(2*np.pi*10*t) + 0.5*np.sin(2*np.pi*50*t) # ノイズ追加 noisy_signal = signal + np.random.normal(0, 0.2, len(t)) # FFTの実行 fft_result = np.fft.fft(noisy_signal) # ローパスフィルタの適用 cutoff_freq = 30 mask = np.abs(freq) < cutoff_freq filtered_fft = fft_result * mask # 逆FFTで信号を復元 filtered_signal = np.real(np.fft.ifft(filtered_fft)) return t, signal, noisy_signal, filtered_signal # 信号処理の実行と結果の表示 t, orig, noisy, filtered = signal_processing_example() print("信号処理完了")
1.4.5 補足：数値計算における複素数の扱い
[編集]科学技術計算では、複素数を扱う際に以下の点に注意が必要です：
- 丸め誤差の影響
- 分岐切断の扱い
- 数値安定性の確保
def numerical_stability_example(): """ 複素数計算における数値安定性の例 """ # 不安定な計算の例 z = 1e15 + 1j w = z + 1 print(f"z = {z}") print(f"w-z = {w-z} # 理論値は1.0") # 改善された計算 z_real = z.real z_imag = z.imag w_improved = complex(z_real + 1, z_imag) print(f"改善された計算: w-z = {w_improved-z}") numerical_stability_example()
1.5 確率・統計
[編集]確率・統計は、データ分析、機械学習、実験データの解析など、幅広い応用場面で必要となる基礎知識です。
1.5.1 確率の基礎
[編集]- 確率空間と確率変数
確率空間は、標本空間Ω、事象の集合F、確率測度Pの三つ組(Ω, F, P)で定義されます。
import numpy as np from scipy import stats # サイコロを振る実験のシミュレーション def dice_experiment(n_trials=10000): """ 公平なサイコロを振る実験のシミュレーション """ results = np.random.randint(1, 7, size=n_trials) probabilities = np.bincount(results)[1:] / n_trials print("実験的確率:") for face, prob in enumerate(probabilities, 1): print(f"目{face}: {prob:.4f}") # カイ二乗検定で一様性を確認 chi2, p_value = stats.chisquare(probabilities) print(f"\nカイ二乗検定:") print(f"統計量: {chi2:.4f}") print(f"p値: {p_value:.4f}") dice_experiment()
1.5.2 確率分布と統計量
[編集]- 代表的な確率分布
def distribution_examples(): """ 代表的な確率分布のサンプリングと可視化 """ np.random.seed(42) n_samples = 1000 # 正規分布 normal_samples = np.random.normal(loc=0, scale=1, size=n_samples) # 指数分布 exponential_samples = np.random.exponential(scale=1, size=n_samples) # 二項分布 binomial_samples = np.random.binomial(n=10, p=0.3, size=n_samples) # 基本統計量の計算 def compute_statistics(samples, name): print(f"\n{name}の統計量:") print(f"平均: {np.mean(samples):.4f}") print(f"分散: {np.var(samples):.4f}") print(f"標準偏差: {np.std(samples):.4f}") print(f"中央値: {np.median(samples):.4f}") compute_statistics(normal_samples, "正規分布") compute_statistics(exponential_samples, "指数分布") compute_statistics(binomial_samples, "二項分布") distribution_examples()
1.5.3 統計的推定と検定
[編集]- 区間推定と仮説検定
def statistical_inference_example(): """ 統計的推定と検定の例 """ # データ生成 np.random.seed(42) population_mean = 100 population_std = 15 sample_size = 30 sample = np.random.normal(population_mean, population_std, sample_size) # 平均の信頼区間 confidence_level = 0.95 t_stat = stats.t.ppf((1 + confidence_level) / 2, df=sample_size-1) sample_mean = np.mean(sample) sample_std = np.std(sample, ddof=1) margin_error = t_stat * (sample_std / np.sqrt(sample_size)) ci_lower = sample_mean - margin_error ci_upper = sample_mean + margin_error print(f"{confidence_level*100}%信頼区間:") print(f"[{ci_lower:.2f}, {ci_upper:.2f}]") # 仮説検定 null_mean = 95 # 帰無仮説の平均値 t_stat, p_value = stats.ttest_1samp(sample, null_mean) print(f"\n一標本のt検定:") print(f"帰無仮説: μ = {null_mean}") print(f"t統計量: {t_stat:.4f}") print(f"p値: {p_value:.4f}") statistical_inference_example()
1.5.4 回帰分析
[編集]- 単回帰と重回帰
from sklearn.linear_model import LinearRegression from sklearn.model_selection import train_test_split from sklearn.metrics import r2_score, mean_squared_error def regression_analysis_example(): """ 回帰分析の実装例 """ # データ生成 np.random.seed(42) n_samples = 100 # 説明変数 X = np.random.rand(n_samples, 2) # 2つの特徴量 # 目的変数（ノイズを含む） y = 2*X[:, 0] + 3*X[:, 1] + np.random.normal(0, 0.1, n_samples) # データの分割 X_train, X_test, y_train, y_test = train_test_split( X, y, test_size=0.2, random_state=42 ) # モデルの学習 model = LinearRegression() model.fit(X_train, y_train) # 予測と評価 y_pred = model.predict(X_test) print("回帰係数:") for i, coef in enumerate(model.coef_): print(f"β{i+1}: {coef:.4f}") print(f"切片: {model.intercept_:.4f}") print(f"\n決定係数 R²: {r2_score(y_test, y_pred):.4f}") print(f"平均二乗誤差: {mean_squared_error(y_test, y_pred):.4f}") regression_analysis_example()
1.5.5 応用例：ブートストラップ法による信頼区間推定
[編集]def bootstrap_example(): """ ブートストラップ法による中央値の信頼区間推定 """ # データ生成 np.random.seed(42) data = np.random.lognormal(0, 0.5, size=100) # ブートストラップサンプリング n_bootstrap = 1000 bootstrap_medians = np.zeros(n_bootstrap) for i in range(n_bootstrap): bootstrap_sample = np.random.choice(data, size=len(data), replace=True) bootstrap_medians[i] = np.median(bootstrap_sample) # 信頼区間の計算（95%） ci_lower = np.percentile(bootstrap_medians, 2.5) ci_upper = np.percentile(bootstrap_medians, 97.5) print("中央値のブートストラップ信頼区間:") print(f"点推定値: {np.median(data):.4f}") print(f"95%信頼区間: [{ci_lower:.4f}, {ci_upper:.4f}]") bootstrap_example()
第2章 数値計算の基礎
[編集]2.1 数値表現と誤差
[編集]計算機における数値計算では、有限の精度で実数を表現する必要があり、そこから生じる誤差を適切に理解し管理することが重要です。
2.1.1 浮動小数点数の表現
[編集]IEEE 754規格による浮動小数点数の表現について理解しましょう。
import sys import numpy as np def floating_point_example(): """ 浮動小数点数の特性を確認する例 """ print(f"機械イプシロン: {np.finfo(float).eps}") print(f"最小の正規化数: {np.finfo(float).tiny}") print(f"最大の有限数: {np.finfo(float).max}") # 丸め誤差の例 x = 0.1 + 0.2 print(f"\n0.1 + 0.2 = {x}") print(f"0.3との差: {x - 0.3}") # 非正規化数の例 x = np.float64(1.0) for _ in range(1074): x /= 2 print(f"x = {x:e}", end='\r') print(f"\n最小の非正規化数: {x:e}") floating_point_example()
2.1.2 数値誤差の種類と対策
[編集]def error_analysis_example(): """ 数値計算における主な誤差の例と対策 """ # 桁落ち（相対的に大きな数からの小さな数の減算） def catastrophic_cancellation(): x = 1.0 delta = 1e-16 bad_result = (x + delta) - x print(f"桁落ちの例:") print(f"期待値: {delta}") print(f"実際の結果: {bad_result}") # 対策: 数式の変形 better_result = delta print(f"改善後の結果: {better_result}") # 丸め誤差の蓄積 def error_accumulation(): sum1 = sum([0.1] * 10) sum2 = 0.1 * 10 print(f"\n丸め誤差の蓄積:") print(f"繰り返し加算: {sum1}") print(f"乗算: {sum2}") print(f"差: {sum1 - sum2}") # 対策: Kahan summation algorithm def kahan_sum(numbers): s = 0.0 c = 0.0 # 補償項 for num in numbers: y = num - c t = s + y c = (t - s) - y s = t return s sum3 = kahan_sum([0.1] * 10) print(f"Kahan summation結果: {sum3}") catastrophic_cancellation() error_accumulation() error_analysis_example()
2.1.3 条件数と数値安定性
[編集]def condition_number_example(): """ 条件数と数値安定性の関係を示す例 """ # 行列の条件数と連立方程式の解の精度 def matrix_condition(): # 良条件の行列 A1 = np.array([[2.0, 1.0], [1.0, 2.0]]) # 悪条件の行列 A2 = np.array([[2.0, 1.999], [1.999, 2.0]]) b = np.array([1.0, 1.0]) print("条件数の比較:") print(f"良条件行列の条件数: {np.linalg.cond(A1):.2f}") print(f"悪条件行列の条件数: {np.linalg.cond(A2):.2f}") # 摂動を加えた場合の解の変化 delta = 1e-10 b_perturbed = b + delta x1 = np.linalg.solve(A1, b) x1_perturbed = np.linalg.solve(A1, b_perturbed) x2 = np.linalg.solve(A2, b) x2_perturbed = np.linalg.solve(A2, b_perturbed) print("\n摂動に対する解の変化:") print("良条件行列:") print(f"相対誤差: {np.linalg.norm(x1_perturbed - x1) / np.linalg.norm(x1):.2e}") print("悪条件行列:") print(f"相対誤差: {np.linalg.norm(x2_perturbed - x2) / np.linalg.norm(x2):.2e}") matrix_condition() condition_number_example()
2.1.4 実践的な対策と推奨事項
[編集]def numerical_best_practices(): """ 数値計算における推奨プラクティス """ # 1. スケーリング def scaling_example(): # 悪い例 x = np.array([1e-8, 1e8]) print(f"スケーリング前の数値範囲: {np.min(x)} to {np.max(x)}") # 良い例 x_scaled = x / np.max(np.abs(x)) print(f"スケーリング後の数値範囲: {np.min(x_scaled)} to {np.max(x_scaled)}") # 2. 数値的に安定なアルゴリズムの選択 def stable_algorithm_example(): # 悪い例：不安定な二次方程式の解法 def unstable_quadratic(a, b, c): x1 = (-b + np.sqrt(b**2 - 4*a*c)) / (2*a) x2 = (-b - np.sqrt(b**2 - 4*a*c)) / (2*a) return x1, x2 # 良い例：数値的に安定な解法 def stable_quadratic(a, b, c): if b >= 0: x1 = (-b - np.sqrt(b**2 - 4*a*c)) / (2*a) x2 = c / (a*x1) else: x1 = (-b + np.sqrt(b**2 - 4*a*c)) / (2*a) x2 = c / (a*x1) return x1, x2 # 比較 a, b, c = 1.0, 1e10, 1.0 print("\n二次方程式の解法比較:") print(f"不安定な方法: {unstable_quadratic(a, b, c)}") print(f"安定な方法: {stable_quadratic(a, b, c)}") scaling_example() stable_algorithm_example() numerical_best_practices()
2.2 方程式の解法
[編集]2.2.1 非線形方程式の数値解法
[編集]非線形方程式 f(x) = 0 の解を求めるための代表的な手法を実装します。
import numpy as np from typing import Callable, Tuple def newton_method(f: Callable[[float], float], df: Callable[[float], float], x0: float, tol: float = 1e-6, max_iter: int = 100) -> Tuple[float, int]: """ ニュートン法による非線形方程式の求解 Parameters: f: 対象となる関数 df: fの導関数 x0: 初期値 tol: 収束判定の閾値 max_iter: 最大反復回数 Returns: (解, 反復回数) """ x = x0 for i in range(max_iter): fx = f(x) if abs(fx) < tol: return x, i dfx = df(x) if abs(dfx) < 1e-10: # 導関数がゼロに近い場合 raise ValueError("導関数がゼロに近いため計算を中断します") x_new = x - fx/dfx if abs(x_new - x) < tol: return x_new, i + 1 x = x_new raise RuntimeError("最大反復回数に達しました") def secant_method(f: Callable[[float], float], x0: float, x1: float, tol: float = 1e-6, max_iter: int = 100) -> Tuple[float, int]: """ 割線法による非線形方程式の求解 """ for i in range(max_iter): f0, f1 = f(x0), f(x1) if abs(f1) < tol: return x1, i if abs(f1 - f0) < 1e-10: raise ValueError("割線の傾きがゼロに近いため計算を中断します") x_new = x1 - f1 * (x1 - x0)/(f1 - f0) if abs(x_new - x1) < tol: return x_new, i + 1 x0, x1 = x1, x_new raise RuntimeError("最大反復回数に達しました") # 方程式を解く例 def equation_solving_example(): """ 非線形方程式を解く例 """ # 例1: x³ - x - 1 = 0 f1 = lambda x: x**3 - x - 1 df1 = lambda x: 3*x**2 - 1 try: solution_newton, iterations = newton_method(f1, df1, 1.0) print("ニュートン法による解:") print(f"x ≈ {solution_newton:.10f}") print(f"反復回数: {iterations}") print(f"残差: {abs(f1(solution_newton)):.2e}") except Exception as e: print(f"エラー: {e}") # 例2: cos(x) = x f2 = lambda x: np.cos(x) - x try: solution_secant, iterations = secant_method(f2, 0.0, 1.0) print("\n割線法による解:") print(f"x ≈ {solution_secant:.10f}") print(f"反復回数: {iterations}") print(f"残差: {abs(f2(solution_secant)):.2e}") except Exception as e: print(f"エラー: {e}") equation_solving_example()
2.2.2 連立方程式の解法
[編集]def linear_system_solvers(): """ 連立一次方程式の各種解法の実装と比較 """ def gauss_seidel(A: np.ndarray, b: np.ndarray, x0: np.ndarray, tol: float = 1e-6, max_iter: int = 1000) -> Tuple[np.ndarray, int]: """ ガウス・ザイデル法による連立一次方程式の求解 """ n = len(b) x = x0.copy() for it in range(max_iter): x_new = x.copy() for i in range(n): s1 = np.dot(A[i, :i], x_new[:i]) s2 = np.dot(A[i, i+1:], x[i+1:]) x_new[i] = (b[i] - s1 - s2) / A[i, i] if np.allclose(x, x_new, rtol=tol): return x_new, it + 1 x = x_new raise RuntimeError("最大反復回数に達しました") # 例題: 3×3の連立方程式 A = np.array([[4.0, -1.0, 0.0], [-1.0, 4.0, -1.0], [0.0, -1.0, 4.0]]) b = np.array([1.0, 5.0, 0.0]) x0 = np.zeros_like(b) # 1. ガウス・ザイデル法 try: solution_gs, iterations = gauss_seidel(A, b, x0) print("ガウス・ザイデル法による解:") print(f"x = {solution_gs}") print(f"反復回数: {iterations}") print(f"残差ノルム: {np.linalg.norm(A @ solution_gs - b):.2e}") except Exception as e: print(f"エラー: {e}") # 2. NumPyの直接解法との比較 solution_np = np.linalg.solve(A, b) print("\nNumPy線形ソルバーによる解:") print(f"x = {solution_np}") print(f"残差ノルム: {np.linalg.norm(A @ solution_np - b):.2e}") # 3. 条件数と精度の関係 cond = np.linalg.cond(A) print(f"\n行列の条件数: {cond:.2e}") print(f"予想される相対誤差の上限: {cond * np.finfo(float).eps:.2e}") linear_system_solvers()
2.2.3 最適化問題の解法
[編集]def optimization_examples(): """ 最適化問題の解法例 """ from scipy.optimize import minimize # 1. 最急降下法の実装 def gradient_descent(f: Callable[[np.ndarray], float], grad_f: Callable[[np.ndarray], np.ndarray], x0: np.ndarray, learning_rate: float = 0.1, tol: float = 1e-6, max_iter: int = 1000) -> Tuple[np.ndarray, float]: """ 最急降下法による最適化 """ x = x0.copy() f_value = f(x) for i in range(max_iter): grad = grad_f(x) if np.linalg.norm(grad) < tol: return x, f_value x_new = x - learning_rate * grad f_value_new = f(x_new) if abs(f_value_new - f_value) < tol: return x_new, f_value_new x = x_new f_value = f_value_new raise RuntimeError("最大反復回数に達しました") # 例題1: ロゼンブロック関数の最小化 def rosenbrock(x): return (1 - x[0])**2 + 100 * (x[1] - x[0]'''2)'''2 def rosenbrock_grad(x): return np.array([ -2*(1 - x[0]) - 400*x[0]*(x[1] - x[0]**2), 200*(x[1] - x[0]**2) ]) x0 = np.array([-1.0, 1.0]) # 最急降下法による解法 try: solution_gd, min_value = gradient_descent( rosenbrock, rosenbrock_grad, x0 ) print("最急降下法による最適化結果:") print(f"x = {solution_gd}") print(f"f(x) = {min_value:.6f}") except Exception as e: print(f"エラー: {e}") # scipyの最適化ソルバーとの比較 result = minimize(rosenbrock, x0, method='BFGS', jac=rosenbrock_grad) print("\nSciPy BFGSによる最適化結果:") print(f"x = {result.x}") print(f"f(x) = {result.fun:.6f}") print(f"収束: {result.success}") print(f"反復回数: {result.nit}") optimization_examples()
2.3 数値微分・積分
[編集]2.3.1 数値微分の手法
[編集]数値微分では、差分近似を用いて導関数を計算します。精度と計算コストのバランスが重要です。
import numpy as np from typing import Callable from dataclasses import dataclass @dataclass class DerivativeResult: """数値微分の結果を格納するクラス""" value: float error_estimate: float step_size: float def numerical_derivative(f: Callable[[float], float], x: float, method: str = 'central') -> DerivativeResult: """ 数値微分の実装 Parameters: f: 微分する関数 x: 微分を計算する点 method: 'forward'（前進差分）, 'central'（中心差分）, または 'richardson'（リチャードソン外挿） """ eps = np.finfo(float).eps h = np.sqrt(eps) * (1 + abs(x)) # 初期ステップサイズ if method == 'forward': derivative = (f(x + h) - f(x)) / h error_est = abs(h) elif method == 'central': derivative = (f(x + h) - f(x - h)) / (2 * h) error_est = abs(h) ** 2 elif method == 'richardson': # リチャードソン外挿 h2 = h / 2 d1 = (f(x + h) - f(x - h)) / (2 * h) d2 = (f(x + h2) - f(x - h2)) / (2 * h2) derivative = (4 * d2 - d1) / 3 # 4次精度の近似 error_est = abs(d2 - d1) / 3 else: raise ValueError("Unknown method") return DerivativeResult(derivative, error_est, h) def higher_derivatives(f: Callable[[float], float], x: float, order: int = 2) -> float: """ 高次導関数の数値計算 """ if order < 1: raise ValueError("Order must be positive") h = np.cbrt(np.finfo(float).eps) * (1 + abs(x)) if order == 1: return numerical_derivative(f, x, 'central').value elif order == 2: # 中心差分による2階導関数 return (f(x + h) - 2*f(x) + f(x - h)) / (h * h) else: # 再帰的に高次導関数を計算 def first_derivative(y): return numerical_derivative(f, y, 'central').value return higher_derivatives(first_derivative, x, order - 1) # 数値微分の例 def derivative_examples(): """数値微分の実践例""" # テスト関数 def f(x): return np.sin(x) * np.exp(-x/2) # 解析的な導関数 def df(x): return np.exp(-x/2) * (np.cos(x) - 0.5*np.sin(x)) x = 1.0 # 各手法での計算結果 methods = ['forward', 'central', 'richardson'] print("数値微分の比較:") for method in methods: result = numerical_derivative(f, x, method) print(f"\n{method}差分法:") print(f"計算値: {result.value:.10f}") print(f"真値との誤差: {abs(result.value - df(x)):.2e}") print(f"推定誤差: {result.error_estimate:.2e}") print(f"使用ステップサイズ: {result.step_size:.2e}") # 高次導関数の計算 print("\n高次導関数:") for order in range(1, 4): deriv = higher_derivatives(f, x, order) print(f"{order}階導関数: {deriv:.10f}") derivative_examples()
2.3.2 数値積分の手法
[編集]@dataclass class IntegrationResult: """数値積分の結果を格納するクラス""" value: float error_estimate: float n_evaluations: int def adaptive_simpson(f: Callable[[float], float], a: float, b: float, tol: float = 1e-8, max_depth: int = 50) -> IntegrationResult: """ 適応的シンプソン法による数値積分 """ def simpson_rule(f, a, b): """基本的なシンプソン則""" c = (a + b) / 2 h = (b - a) / 6 return h * (f(a) + 4*f(c) + f(b)) def _adaptive_simpson_recurse(a, b, fa, fb, fc, integral, tol, depth): """再帰的な適応的シンプソン法の実装""" c = (a + b) / 2 h = (b - a) / 12 d = (a + c) / 2 e = (c + b) / 2 fd = f(d) fe = f(e) left_integral = h * (fa + 4*fd + fc) right_integral = h * (fc + 4*fe + fb) whole_integral = left_integral + right_integral if depth >= max_depth: return whole_integral, abs(whole_integral - integral) if abs(whole_integral - integral) <= 15 * tol: return whole_integral, abs(whole_integral - integral) / 15 left_result, left_error = _adaptive_simpson_recurse( a, c, fa, fc, fd, left_integral, tol/2, depth + 1) right_result, right_error = _adaptive_simpson_recurse( c, b, fc, fb, fe, right_integral, tol/2, depth + 1) return left_result + right_result, left_error + right_error # 初期評価 c = (a + b) / 2 fa, fb, fc = f(a), f(b), f(c) initial_integral = simpson_rule(f, a, b) result, error = _adaptive_simpson_recurse( a, b, fa, fb, fc, initial_integral, tol, 0) return IntegrationResult(result, error, max_depth) def gauss_quadrature(f: Callable[[float], float], a: float, b: float, n: int = 5) -> IntegrationResult: """ ガウス求積法による数値積分 """ # ガウス・ルジャンドル求積の重みと節点（5点法） x5 = np.array([0.0, ±0.5384693101056831, ±0.9061798459386640]) w5 = np.array([0.5688888888888889, 0.4786286704993665, 0.2369268850561891]) # 区間変換 mid = (b + a) / 2 half_length = (b - a) / 2 x = mid + half_length * x5 # 積分計算 result = half_length * np.sum(w5 * f(x)) # 誤差推定（簡易的） error_est = abs(result) * 1e-12 return IntegrationResult(result, error_est, 5) # 数値積分の例 def integration_examples(): """数値積分の実践例""" # テスト関数 def f(x): return np.exp(-x*x) * np.sin(2*x) # 積分区間 a, b = 0, 2 # 適応的シンプソン法 result_simpson = adaptive_simpson(f, a, b) print("適応的シンプソン法:") print(f"積分値: {result_simpson.value:.10f}") print(f"推定誤差: {result_simpson.error_estimate:.2e}") # ガウス求積法 result_gauss = gauss_quadrature(f, a, b) print("\nガウス求積法:") print(f"積分値: {result_gauss.value:.10f}") print(f"推定誤差: {result_gauss.error_estimate:.2e}") # scipy.integrateとの比較 from scipy import integrate result_scipy, error_scipy = integrate.quad(f, a, b) print("\nSciPy quad:") print(f"積分値: {result_scipy:.10f}") print(f"推定誤差: {error_scipy:.2e}") integration_examples()
2.4 線形方程式の解法
[編集]線形方程式の効率的な解法は科学技術計算の基礎となります。ここでは、直接法と反復法の両方について解説します。
2.4.1 直接法
[編集]直接法は、有限回の演算で厳密解（丸め誤差を除く）を得る方法です。
import numpy as np from typing import Tuple def lu_decomposition(A: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: """ LU分解を行う関数 Parameters: A: 分解する行列 Returns: L: 下三角行列 U: 上三角行列 """ n = len(A) L = np.zeros((n, n)) U = np.zeros((n, n)) for i in range(n): # L の対角成分 L[i, i] = 1.0 # U の i 行目 for j in range(i, n): sum_u = sum(L[i, k] * U[k, j] for k in range(i)) U[i, j] = A[i, j] - sum_u # L の i 列目 for j in range(i + 1, n): sum_l = sum(L[j, k] * U[k, i] for k in range(i)) L[j, i] = (A[j, i] - sum_l) / U[i, i] return L, U def solve_lu(L: np.ndarray, U: np.ndarray, b: np.ndarray) -> np.ndarray: """ LU分解を用いて連立方程式を解く Parameters: L: 下三角行列 U: 上三角行列 b: 右辺ベクトル Returns: x: 解ベクトル """ n = len(b) # 前進代入 (Ly = b) y = np.zeros(n) for i in range(n): y[i] = b[i] - np.dot(L[i, :i], y[:i]) # 後退代入 (Ux = y) x = np.zeros(n) for i in range(n-1, -1, -1): x[i] = (y[i] - np.dot(U[i, i+1:], x[i+1:])) / U[i, i] return x
2.4.2 反復法
[編集]反復法は、初期値から出発して徐々に解に近づけていく方法です。大規模な疎行列に対して特に有効です。
def conjugate_gradient(A: np.ndarray, b: np.ndarray, x0: np.ndarray = None, tol: float = 1e-10, max_iter: int = 1000) -> Tuple[np.ndarray, int]: """ 共役勾配法による連立方程式の求解 Parameters: A: 係数行列（対称正定値） b: 右辺ベクトル x0: 初期値（デフォルトはゼロベクトル） tol: 収束判定閾値 max_iter: 最大反復回数 Returns: x: 解ベクトル iter_count: 反復回数 """ n = len(b) if x0 is None: x = np.zeros(n) else: x = x0.copy() r = b - A @ x p = r.copy() r_norm_sq = r @ r for i in range(max_iter): Ap = A @ p alpha = r_norm_sq / (p @ Ap) x += alpha * p r -= alpha * Ap r_norm_sq_new = r @ r if np.sqrt(r_norm_sq_new) < tol: return x, i + 1 beta = r_norm_sq_new / r_norm_sq r_norm_sq = r_norm_sq_new p = r + beta * p raise RuntimeError("最大反復回数に達しました")
2.4.3 前処理付き反復法
[編集]前処理を適用することで、反復法の収束性を大幅に改善できます。
def incomplete_cholesky(A: np.ndarray) -> np.ndarray: """ 不完全コレスキー分解による前処理行列の生成 Parameters: A: 対称正定値行列 Returns: L: 下三角行列 """ n = len(A) L = np.zeros((n, n)) for i in range(n): for j in range(i + 1): s = A[i, j] for k in range(j): s -= L[i, k] * L[j, k] if i == j: if s <= 0: raise ValueError("行列が正定値でない可能性があります") L[i, j] = np.sqrt(s) else: L[i, j] = s / L[j, j] return L def preconditioned_cg(A: np.ndarray, b: np.ndarray, M: np.ndarray, x0: np.ndarray = None, tol: float = 1e-10, max_iter: int = 1000) -> Tuple[np.ndarray, int]: """ 前処理付き共役勾配法 Parameters: A: 係数行列 b: 右辺ベクトル M: 前処理行列 x0: 初期値 tol: 収束判定閾値 max_iter: 最大反復回数 Returns: x: 解ベクトル iter_count: 反復回数 """ n = len(b) if x0 is None: x = np.zeros(n) else: x = x0.copy() r = b - A @ x z = np.linalg.solve(M, r) p = z.copy() rz = r @ z for i in range(max_iter): Ap = A @ p alpha = rz / (p @ Ap) x += alpha * p r -= alpha * Ap if np.linalg.norm(r) < tol: return x, i + 1 z = np.linalg.solve(M, r) rz_new = r @ z beta = rz_new / rz rz = rz_new p = z + beta * p raise RuntimeError("最大反復回数に達しました")
2.5 固有値計算
[編集]固有値問題は、振動解析、安定性解析、主成分分析など、様々な応用で重要な役割を果たします。
2.5.1 べき乗法
[編集]べき乗法は最大固有値とそれに対応する固有ベクトルを求める最も基本的な方法です。
def power_method(A: np.ndarray, tol: float = 1e-10, max_iter: int = 1000) -> Tuple[float, np.ndarray]: """ べき乗法による最大固有値・固有ベクトルの計算 Parameters: A: 正方行列 tol: 収束判定閾値 max_iter: 最大反復回数 Returns: eigenvalue: 最大固有値 eigenvector: 対応する固有ベクトル """ n = len(A) x = np.random.rand(n) x = x / np.linalg.norm(x) lambda_old = 0 for i in range(max_iter): # 行列とベクトルの積 y = A @ x # レイリー商による固有値の計算 lambda_new = x @ y # ベクトルの正規化 x = y / np.linalg.norm(y) # 収束判定 if abs(lambda_new - lambda_old) < tol: return lambda_new, x lambda_old = lambda_new raise RuntimeError("最大反復回数に達しました")
第3章 信号処理
[編集]3.1 フーリエ解析
[編集]フーリエ解析は信号処理の基礎となる重要な手法です。時間領域の信号を周波数領域で解析することで、信号の特性を理解し、効果的な処理が可能になります。
3.1.1 離散フーリエ変換（DFT）の基礎
[編集]import numpy as np from typing import Tuple import matplotlib.pyplot as plt def dft(x: np.ndarray) -> np.ndarray: """ 離散フーリエ変換の直接計算による実装 Parameters: x: 入力信号 Returns: X: 周波数領域の信号 """ N = len(x) X = np.zeros(N, dtype=complex) for k in range(N): for n in range(N): X[k] += x[n] * np.exp(-2j * np.pi * k * n / N) return X def idft(X: np.ndarray) -> np.ndarray: """ 逆離散フーリエ変換の直接計算による実装 Parameters: X: 周波数領域の信号 Returns: x: 時間領域の信号 """ N = len(X) x = np.zeros(N, dtype=complex) for n in range(N): for k in range(N): x[n] += X[k] * np.exp(2j * np.pi * k * n / N) return x / N def analyze_signal(t: np.ndarray, x: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: """ 信号のスペクトル解析を行う Parameters: t: 時間配列 x: 信号配列 Returns: freq: 周波数配列 spectrum: パワースペクトル """ N = len(t) dt = t[1] - t[0] # サンプリング間隔 # FFTの計算 X = np.fft.fft(x) # 周波数軸の生成 freq = np.fft.fftfreq(N, dt) # パワースペクトルの計算 spectrum = np.abs(X) ** 2 return freq, spectrum
3.1.2 窓関数と周波数漏れ
[編集]def apply_window(x: np.ndarray, window_type: str = 'hanning') -> np.ndarray: """ 信号に窓関数を適用する Parameters: x: 入力信号 window_type: 窓関数の種類 ('hanning', 'hamming', 'blackman') Returns: x_windowed: 窓関数適用後の信号 """ N = len(x) if window_type == 'hanning': window = 0.5 * (1 - np.cos(2 * np.pi * np.arange(N) / (N - 1))) elif window_type == 'hamming': window = 0.54 - 0.46 * np.cos(2 * np.pi * np.arange(N) / (N - 1)) elif window_type == 'blackman': window = (0.42 - 0.5 * np.cos(2 * np.pi * np.arange(N) / (N - 1)) + 0.08 * np.cos(4 * np.pi * np.arange(N) / (N - 1))) else: raise ValueError("未対応の窓関数です") return x * window def spectral_analysis_with_window(t: np.ndarray, x: np.ndarray, window_type: str = 'hanning') -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """ 窓関数を適用してスペクトル解析を行う Parameters: t: 時間配列 x: 信号配列 window_type: 窓関数の種類 Returns: freq: 周波数配列 spectrum_raw: 生のパワースペクトル spectrum_windowed: 窓関数適用後のパワースペクトル """ # 生信号のスペクトル解析 freq, spectrum_raw = analyze_signal(t, x) # 窓関数の適用 x_windowed = apply_window(x, window_type) # 窓関数適用後の信号のスペクトル解析 _, spectrum_windowed = analyze_signal(t, x_windowed) return freq, spectrum_raw, spectrum_windowed
3.1.3 短時間フーリエ変換（STFT）
[編集]def stft(x: np.ndarray, window_size: int, hop_size: int, window_type: str = 'hanning') -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """ 短時間フーリエ変換の実装 Parameters: x: 入力信号 window_size: 窓サイズ hop_size: フレームシフト量 window_type: 窓関数の種類 Returns: times: 時間フレーム frequencies: 周波数ビン spectrogram: スペクトログラム """ # 窓関数の生成 window = apply_window(np.ones(window_size), window_type) # フレーム数の計算 n_frames = 1 + (len(x) - window_size) // hop_size # スペクトログラムの初期化 spectrogram = np.zeros((window_size // 2 + 1, n_frames), dtype=complex) # STFTの計算 for i in range(n_frames): start = i * hop_size frame = x[start:start + window_size] * window spectrum = np.fft.rfft(frame) spectrogram[:, i] = spectrum # 時間軸と周波数軸の生成 times = np.arange(n_frames) * hop_size frequencies = np.fft.rfftfreq(window_size) return times, frequencies, spectrogram def istft(spectrogram: np.ndarray, window_size: int, hop_size: int, window_type: str = 'hanning') -> np.ndarray: """ 逆短時間フーリエ変換の実装 Parameters: spectrogram: スペクトログラム window_size: 窓サイズ hop_size: フレームシフト量 window_type: 窓関数の種類 Returns: x_reconstructed: 再構成された信号 """ # 窓関数の生成 window = apply_window(np.ones(window_size), window_type) n_frames = spectrogram.shape[1] expected_signal_len = hop_size * (n_frames - 1) + window_size # 信号の再構成 x_reconstructed = np.zeros(expected_signal_len) window_sum = np.zeros(expected_signal_len) for i in range(n_frames): start = i * hop_size frame = np.fft.irfft(spectrogram[:, i]) x_reconstructed[start:start + window_size] += frame * window window_sum[start:start + window_size] += window ** 2 # オーバーラップアド処理 non_zero = window_sum > 1e-10 x_reconstructed[non_zero] /= window_sum[non_zero] return x_reconstructed
3.2 ディジタルフィルタ
[編集]3.2.1 FIRフィルタの設計
[編集]def design_fir_filter(cutoff: float, filter_order: int, filter_type: str = 'lowpass') -> np.ndarray: """ 窓関数法によるFIRフィルタの設計 Parameters: cutoff: カットオフ周波数 (0 < cutoff < 0.5) filter_order: フィルタの次数 filter_type: フィルタの種類 ('lowpass' or 'highpass') Returns: h: フィルタ係数 """ if not 0 < cutoff < 0.5: raise ValueError("カットオフ周波数は0と0.5の間である必要があります") # インパルス応答の計算 n = np.arange(filter_order + 1) omega_c = 2 * np.pi * cutoff if filter_type == 'lowpass': h = np.sinc(omega_c * (n - filter_order/2) / np.pi) elif filter_type == 'highpass': h = np.sinc(n - filter_order/2) - np.sinc(omega_c * (n - filter_order/2) / np.pi) else: raise ValueError("未対応のフィルタ種類です") # ハミング窓の適用 window = 0.54 - 0.46 * np.cos(2 * np.pi * n / filter_order) h = h * window # 正規化 h = h / np.sum(h) return h def apply_fir_filter(x: np.ndarray, h: np.ndarray) -> np.ndarray: """ FIRフィルタの適用 Parameters: x: 入力信号 h: フィルタ係数 Returns: y: フィルタリング後の信号 """ return np.convolve(x, h, mode='same')
3.2.2 IIRフィルタの設計
[編集]from scipy import signal def design_iir_butterworth(cutoff: float, order: int, filter_type: str = 'lowpass') -> Tuple[np.ndarray, np.ndarray]: """ バターワースIIRフィルタの設計 Parameters: cutoff: カットオフ周波数 (0 < cutoff < 0.5) order: フィルタの次数 filter_type: フィルタの種類 Returns: b: 分子の係数 a: 分母の係数 """ nyquist = 0.5 cutoff_normalized = cutoff / nyquist return signal.butter(order, cutoff_normalized, btype=filter_type) def apply_iir_filter(x: np.ndarray, b: np.ndarray, a: np.ndarray, initial_conditions: np.ndarray = None) -> np.ndarray: """ IIRフィルタの適用 Parameters: x: 入力信号 b: 分子の係数 a: 分母の係数 initial_conditions: 初期状態 Returns: y: フィルタリング後の信号 """ return signal.lfilter(b, a, x, zi=initial_conditions)
3.3 画像処理の基礎
[編集]3.3.1 空間フィルタリング
[編集]def spatial_filter(image: np.ndarray, kernel: np.ndarray) -> np.ndarray: """ 2次元畳み込みによる空間フィルタリング Parameters: image: 入力画像 kernel: フィルタカーネル Returns: filtered_image: フィルタリング後の画像 """ # パディングサイズの計算 pad_h = kernel.shape[0] // 2 pad_w = kernel.shape[1] // 2 # パディング処理 padded_image = np.pad(image, ((pad_h, pad_h), (pad_w, pad_w)), mode='reflect') # 出力画像の初期化 filtered_image = np.zeros_like(image, dtype=float) # 畳み込み演算 for i in range(image.shape[0]): for j in range(image.shape[1]): window = padded_image[i:i + kernel.shape[0], j:j + kernel.shape[1]] filtered_image[i, j] = np.sum(window * kernel) return filtered_image def create_gaussian_kernel(size: int, sigma: float) -> np.ndarray: """ ガウシアンフィルタカーネルの生成 Parameters: size: カーネルサイズ sigma: 標準偏差 Returns: kernel: ガウシアンカーネル """ x, y = np.meshgrid(np.linspace(-size//2, size//2, size), np.linspace(-size//2, size//2, size)) kernel = np.exp(-(x'''2 + y'''2) / (2 * sigma**2)) return kernel / kernel.sum()
3.3.2 エッジ検出
[編集]import numpy as np from typing import Tuple def sobel_edge_detection(image: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """ Sobelオペレータによるエッジ検出 Parameters: image: 入力画像（グレースケール） Returns: magnitude: エッジの強度 direction: エッジの方向 edges: 二値化されたエッジ画像 """ # Sobelカーネルの定義 kernel_x = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]]) kernel_y = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]]) # 勾配の計算 gradient_x = spatial_filter(image, kernel_x) gradient_y = spatial_filter(image, kernel_y) # エッジの強度と方向の計算 magnitude = np.sqrt(gradient_x'''2 + gradient_y'''2) direction = np.arctan2(gradient_y, gradient_x) # 閾値処理によるエッジの二値化 threshold = 0.5 * np.max(magnitude) edges = magnitude > threshold return magnitude, direction, edges def canny_edge_detection(image: np.ndarray, low_threshold: float = 0.1, high_threshold: float = 0.3) -> np.ndarray: """ Cannyエッジ検出器の実装 Parameters: image: 入力画像 low_threshold: 低閾値（最大値に対する比率） high_threshold: 高閾値（最大値に対する比率） Returns: edges: 検出されたエッジ """ # ステップ1: ガウシアンフィルタによるノイズ除去 smoothed = spatial_filter(image, create_gaussian_kernel(5, 1.4)) # ステップ2: 勾配の計算 magnitude, direction, _ = sobel_edge_detection(smoothed) # ステップ3: 非最大値抑制 suppressed = non_maximum_suppression(magnitude, direction) # ステップ4: ヒステリシス閾値処理 high = high_threshold * np.max(suppressed) low = low_threshold * np.max(suppressed) return hysteresis_thresholding(suppressed, low, high) def non_maximum_suppression(magnitude: np.ndarray, direction: np.ndarray) -> np.ndarray: """ 非最大値抑制の実装 """ height, width = magnitude.shape suppressed = np.zeros_like(magnitude) # 勾配方向を0°、45°、90°、135°に量子化 direction = np.rad2deg(direction) % 180 for i in range(1, height - 1): for j in range(1, width - 1): if direction[i, j] < 22.5 or direction[i, j] >= 157.5: neighbors = [magnitude[i, j-1], magnitude[i, j+1]] elif 22.5 <= direction[i, j] < 67.5: neighbors = [magnitude[i-1, j-1], magnitude[i+1, j+1]] elif 67.5 <= direction[i, j] < 112.5: neighbors = [magnitude[i-1, j], magnitude[i+1, j]] else: neighbors = [magnitude[i-1, j+1], magnitude[i+1, j-1]] if magnitude[i, j] >= max(neighbors): suppressed[i, j] = magnitude[i, j] return suppressed
3.3.3 特徴点検出
[編集]def harris_corner_detector(image: np.ndarray, k: float = 0.04, threshold: float = 0.01) -> np.ndarray: """ Harrisコーナー検出器の実装 Parameters: image: 入力画像 k: Harrisレスポンス関数のパラメータ threshold: コーナー検出の閾値 Returns: corners: コーナー点の位置（True/False） """ # 勾配の計算 dx = spatial_filter(image, np.array([[-1, 0, 1]])) dy = spatial_filter(image, np.array([[-1, 0, 1]]).T) # 構造テンソルの要素 Ixx = spatial_filter(dx**2, create_gaussian_kernel(3, 1.0)) Iyy = spatial_filter(dy**2, create_gaussian_kernel(3, 1.0)) Ixy = spatial_filter(dx*dy, create_gaussian_kernel(3, 1.0)) # Harrisレスポンスの計算 det = Ixx * Iyy - Ixy**2 trace = Ixx + Iyy response = det - k * trace**2 # 閾値処理と非最大値抑制 corners = np.zeros_like(image, dtype=bool) response_max = np.max(response) threshold_abs = threshold * response_max for i in range(1, image.shape[0]-1): for j in range(1, image.shape[1]-1): if response[i,j] > threshold_abs: window = response[i-1:i+2, j-1:j+2] if response[i,j] == np.max(window): corners[i,j] = True return corners
3.4 時系列解析
[編集]3.4.1 自己相関分析
[編集]def autocorrelation(x: np.ndarray, max_lag: int = None) -> np.ndarray: """ 時系列データの自己相関関数を計算 Parameters: x: 入力時系列データ max_lag: 最大ラグ（デフォルトはデータ長の半分） Returns: acf: 自己相関関数 """ if max_lag is None: max_lag = len(x) // 2 # 平均を引く x_centered = x - np.mean(x) # 自己相関の計算 acf = np.zeros(max_lag + 1) variance = np.var(x_centered) for lag in range(max_lag + 1): acf[lag] = np.sum(x_centered[lag:] * x_centered[:-lag if lag > 0 else None]) # 正規化 return acf / (len(x) * variance) def partial_autocorrelation(x: np.ndarray, max_lag: int = None) -> np.ndarray: """ 偏自己相関関数の計算（Durbin-Levinsonアルゴリズム） Parameters: x: 入力時系列データ max_lag: 最大ラグ Returns: pacf: 偏自己相関関数 """ if max_lag is None: max_lag = len(x) // 2 # 自己相関関数の計算 acf = autocorrelation(x, max_lag) # 偏自己相関係数の初期化 pacf = np.zeros(max_lag + 1) pacf[0] = 1.0 # Durbin-Levinsonアルゴリズム for k in range(1, max_lag + 1): phi = np.zeros(k) phi[k-1] = acf[k] for j in range(k-1): phi[j] = pacf[j] for j in range(k-1): phi[j] -= pacf[k-1] * pacf[k-2-j] pacf[k-1] = phi[k-1] return pacf
3.4.2 スペクトル推定
[編集]def welch_periodogram(x: np.ndarray, window_size: int = None, overlap: float = 0.5, window_type: str = 'hanning') -> Tuple[np.ndarray, np.ndarray]: """ Welch法によるパワースペクトル密度の推定 Parameters: x: 入力時系列データ window_size: 窓サイズ overlap: オーバーラップ率（0～1） window_type: 窓関数の種類 Returns: frequencies: 周波数配列 psd: パワースペクトル密度 """ if window_size is None: window_size = len(x) // 8 # オーバーラップするサンプル数 hop_size = int(window_size * (1 - overlap)) # 窓関数の適用 window = apply_window(np.ones(window_size), window_type) # セグメントの数の計算 n_segments = (len(x) - window_size) // hop_size + 1 # スペクトルの初期化 spectrum = np.zeros(window_size // 2 + 1) # 各セグメントのスペクトルの計算と平均化 for i in range(n_segments): start = i * hop_size segment = x[start:start + window_size] * window segment_spectrum = np.abs(np.fft.rfft(segment))**2 spectrum += segment_spectrum # 正規化 spectrum /= n_segments # 周波数軸の生成 frequencies = np.fft.rfftfreq(window_size) return frequencies, spectrum
3.5 ウェーブレット変換
[編集]3.5.1 連続ウェーブレット変換
[編集]def morlet_wavelet(t: np.ndarray, omega0: float = 6.0) -> np.ndarray: """ モルレーウェーブレットの生成 Parameters: t: 時間配列 omega0: 中心周波数 Returns: psi: ウェーブレット関数 """ return np.pi**(-0.25) * np.exp(1j * omega0 * t) * np.exp(-t**2 / 2) def cwt(x: np.ndarray, scales: np.ndarray, dt: float = 1.0) -> np.ndarray: """ 連続ウェーブレット変換の実装 Parameters: x: 入力信号 scales: スケールパラメータの配列 dt: サンプリング間隔 Returns: W: ウェーブレット係数（時間-スケール表現） """ n = len(x) W = np.zeros((len(scales), n), dtype=complex) # 各スケールでのウェーブレット変換 for i, scale in enumerate(scales): # ウェーブレットの生成 t = np.arange((-n//2 + 1), (n//2)) * dt wavelet = morlet_wavelet(t/scale) # 畳み込み（周波数領域で計算） x_pad = np.pad(x, (n//2, n//2), mode='reflect') convolution = np.fft.ifft(np.fft.fft(x_pad) * np.fft.fft(wavelet)) W[i, :] = convolution[n//2:-n//2] return W
第4章 シミュレーション技法
[編集]4.1 常微分方程式の解法
[編集]常微分方程式は物理現象や工学システムのモデリングにおいて重要な役割を果たします。ここでは数値解法の基本的な手法を実装します。
4.1.1 オイラー法とルンゲ・クッタ法
[編集]import numpy as np from typing import Tuple, Callable def euler_method(f: Callable, y0: float, t_span: Tuple[float, float], h: float) -> Tuple[np.ndarray, np.ndarray]: """ オイラー法による常微分方程式の数値解法 Parameters: f: dy/dt = f(t, y) の右辺を表す関数 y0: 初期値 t_span: 時間の範囲 (t_start, t_end) h: 時間ステップ Returns: t: 時間配列 y: 解の配列 """ t_start, t_end = t_span n_steps = int((t_end - t_start) / h) t = np.linspace(t_start, t_end, n_steps + 1) y = np.zeros(n_steps + 1) y[0] = y0 for i in range(n_steps): y[i + 1] = y[i] + h * f(t[i], y[i]) return t, y def runge_kutta4(f: Callable, y0: float, t_span: Tuple[float, float], h: float) -> Tuple[np.ndarray, np.ndarray]: """ 4次のルンゲ・クッタ法による常微分方程式の数値解法 Parameters: f: dy/dt = f(t, y) の右辺を表す関数 y0: 初期値 t_span: 時間の範囲 (t_start, t_end) h: 時間ステップ Returns: t: 時間配列 y: 解の配列 """ t_start, t_end = t_span n_steps = int((t_end - t_start) / h) t = np.linspace(t_start, t_end, n_steps + 1) y = np.zeros(n_steps + 1) y[0] = y0 for i in range(n_steps): k1 = f(t[i], y[i]) k2 = f(t[i] + h/2, y[i] + h*k1/2) k3 = f(t[i] + h/2, y[i] + h*k2/2) k4 = f(t[i] + h, y[i] + h*k3) y[i + 1] = y[i] + h * (k1 + 2*k2 + 2*k3 + k4) / 6 return t, y
4.1.2 適応的ステップサイズ制御
[編集]def adaptive_rk4(f: Callable, y0: float, t_span: Tuple[float, float], h0: float, tol: float = 1e-6) -> Tuple[np.ndarray, np.ndarray]: """ 適応的なステップサイズ制御を行う4次のルンゲ・クッタ法 Parameters: f: dy/dt = f(t, y) の右辺を表す関数 y0: 初期値 t_span: 時間の範囲 (t_start, t_end) h0: 初期時間ステップ tol: 許容誤差 Returns: t: 時間配列 y: 解の配列 """ t_start, t_end = t_span t = [t_start] y = [y0] h = h0 while t[-1] < t_end: # 2つの異なるステップサイズでの解を計算 k1 = f(t[-1], y[-1]) k2 = f(t[-1] + h/2, y[-1] + h*k1/2) k3 = f(t[-1] + h/2, y[-1] + h*k2/2) k4 = f(t[-1] + h, y[-1] + h*k3) y1 = y[-1] + h * (k1 + 2*k2 + 2*k3 + k4) / 6 # 半分のステップサイズでの解 h_half = h/2 k1 = f(t[-1], y[-1]) k2 = f(t[-1] + h_half/2, y[-1] + h_half*k1/2) k3 = f(t[-1] + h_half/2, y[-1] + h_half*k2/2) k4 = f(t[-1] + h_half, y[-1] + h_half*k3) y_half = y[-1] + h_half * (k1 + 2*k2 + 2*k3 + k4) / 6 k1 = f(t[-1] + h_half, y_half) k2 = f(t[-1] + h_half + h_half/2, y_half + h_half*k1/2) k3 = f(t[-1] + h_half + h_half/2, y_half + h_half*k2/2) k4 = f(t[-1] + h, y_half + h_half*k3) y2 = y_half + h_half * (k1 + 2*k2 + 2*k3 + k4) / 6 # 誤差の推定 error = abs(y2 - y1) # ステップサイズの調整 if error < tol: t.append(t[-1] + h) y.append(y2) h = min(h * (tol/error)**(1/4), 2*h) # ステップサイズの増加を制限 else: h = h * (tol/error)**(1/4) return np.array(t), np.array(y)
4.2 偏微分方程式の解法
[編集]偏微分方程式は、熱伝導や波動現象など、時間と空間の両方に依存する現象を記述するのに使用されます。ここでは、差分法による基本的な解法を実装します。
4.2.1 熱伝導方程式
[編集]def heat_equation_1d(u0: np.ndarray, dx: float, dt: float, D: float, n_steps: int) -> np.ndarray: """ 1次元熱伝導方程式の陽的差分法による解法 ∂u/∂t = D ∂²u/∂x² Parameters: u0: 初期温度分布 dx: 空間刻み幅 dt: 時間刻み幅 D: 熱拡散係数 n_steps: 時間ステップ数 Returns: u: 各時刻での温度分布 """ nx = len(u0) u = np.zeros((n_steps + 1, nx)) u[0] = u0 # 安定性条件のチェック if D * dt / (dx**2) > 0.5: raise ValueError("数値不安定性の可能性があります") for n in range(n_steps): for i in range(1, nx-1): u[n+1,i] = u[n,i] + D * dt / (dx**2) * (u[n,i+1] - 2*u[n,i] + u[n,i-1]) # 境界条件（固定端） u[n+1,0] = u[n+1,1] u[n+1,-1] = u[n+1,-2] return u
4.2.2 波動方程式
[編集]def wave_equation_1d(u0: np.ndarray, v0: np.ndarray, dx: float, dt: float, c: float, n_steps: int) -> np.ndarray: """ 1次元波動方程式の差分法による解法 ∂²u/∂t² = c² ∂²u/∂x² Parameters: u0: 初期変位 v0: 初期速度 dx: 空間刻み幅 dt: 時間刻み幅 c: 波速 n_steps: 時間ステップ数 Returns: u: 各時刻での変位分布 """ nx = len(u0) u = np.zeros((n_steps + 1, nx)) u[0] = u0 # 安定性条件のチェック（CFL条件） if c * dt / dx > 1: raise ValueError("CFL条件を満たしていません") # 1ステップ目の計算 for i in range(1, nx-1): u[1,i] = u0[i] + v0[i] * dt + (c*dt/dx)**2 * (u0[i+1] - 2*u0[i] + u0[i-1]) / 2 # 境界条件 u[1,0] = 0 u[1,-1] = 0 # 2ステップ目以降の計算 for n in range(1, n_steps): for i in range(1, nx-1): u[n+1,i] = 2*u[n,i] - u[n-1,i] + (c*dt/dx)**2 * (u[n,i+1] - 2*u[n,i] + u[n,i-1]) # 境界条件（固定端） u[n+1,0] = 0 u[n+1,-1] = 0 return u
4.3 モンテカルロ法
[編集]モンテカルロ法は、確率的なサンプリングを用いて数値計算や最適化を行う手法です。高次元の問題や解析的に解くことが困難な問題に特に有効です。
4.3.1 モンテカルロ積分
[編集]def monte_carlo_integration(f: Callable, a: float, b: float, n_samples: int = 10000) -> Tuple[float, float]: """ モンテカルロ積分 Parameters: f: 被積分関数 a: 積分区間の下限 b: 積分区間の上限 n_samples: サンプル数 Returns: integral: 積分値の推定値 error: 標準誤差 """ x = np.random.uniform(a, b, n_samples) y = f(x) integral = (b - a) * np.mean(y) error = (b - a) * np.std(y) / np.sqrt(n_samples) return integral, error
4.3.2 マルコフ連鎖モンテカルロ法
[編集]def metropolis_hastings(p: Callable, q: Callable, q_sample: Callable, n_samples: int) -> np.ndarray: """ メトロポリス・ヘイスティングス法によるMCMCサンプリング Parameters: p: 目的の確率分布 q: 提案分布の確率密度関数 q_sample: 提案分布からのサンプリング関数 n_samples: 生成するサンプル数 Returns: samples: 生成されたサンプル """ samples = np.zeros(n_samples) current = q_sample() # 初期値 for i in range(n_samples): # 新しい状態の提案 proposal = q_sample() # 受理確率の計算 acceptance_ratio = min(1, (p(proposal) * q(current, proposal)) / (p(current) * q(proposal, current))) # 受理判定 if np.random.random() < acceptance_ratio: current = proposal samples[i] = current return samples