科学技術計算.jl

はじめに

本書は、科学技術計算に携わる研究者、技術者、学生のための実践的なハンドブックです。基礎的な数学から最新の機械学習手法まで、幅広い内容を網羅しています。各章には具体的な計算例とJuliaによるプログラミング例を収録し、実務での即戦力となる知識の習得を目指しています。

第1章基礎数学

1.1 代数学の基礎

数学的な記述の基礎となる代数学について解説します。特に、科学技術計算で頻繁に現れる多項式、行列、ベクトルの取り扱いに重点を置きます。

1.1.1 集合と写像

集合論の基本概念から始めましょう。集合とは、ものの集まりを表す最も基本的な概念です。

集合: 集合とは、明確な基準によって、ものの属する・属さないが判定できる、ものの集まりです。

例えば、以下のような集合が考えられます：

自然数の集合 $\mathbb {N} =\{1,2,3,\dots \}$
整数の集合 $\mathbb {Z} =\{\dots ,-2,-1,0,1,2,\dots \}$
実数の集合 $\mathbb {R}$

集合の演算

集合 $A$ , $B$ に対して、以下の基本的な演算が定義されます：

和集合 $A\cup B$ ： $A$ または $B$ のいずれかに属する要素の集合
積集合 $A\cap B$ ： $A$ と $B$ の両方に属する要素の集合
差集合 $A\setminus B$ ： $A$ に属し、 $B$ に属さない要素の集合

1.1.2 群・環・体

代数的構造の基本概念である群、環、体について説明します。

群: 集合 $G$ と二項演算 $\cdot$ が与えられているとき、以下の4つの条件を満たすとき、 $(G,\cdot )$ を群といいます：

結合法則：任意の $a,b,c\in G$ に対して $(a\cdot b)\cdot c=a\cdot (b\cdot c)$
単位元の存在： $e\in G$ が存在して、任意の $a\in G$ に対して $a\cdot e=e\cdot a=a$
逆元の存在：任意の $a\in G$ に対して、 $a\cdot b=b\cdot a=e$ となる $b\in G$ が存在
閉じている：任意の $a,b\in G$ に対して $a\cdot b\in G$

実践例：

using LinearAlgebra

# 回転行列の群の例
rotation_matrix(theta) =[cos(theta) -sin(theta); sin(theta) cos(theta)]

# 2つの回転の合成
theta1 = π/4  # 45度
theta2 = π/3  # 60度
R1 = rotation_matrix(theta1)
R2 = rotation_matrix(theta2)

# 結合法則の確認
R3 = rotation_matrix(theta1 + theta2)
println(isapprox(R1 * R2, R3))  # true

1.2 微分積分学

1.2.1 微分の基礎

微分は、関数の局所的な変化率を表す重要な概念です。科学技術計算において、物理現象の記述や最適化問題の解法に不可欠です。

導関数の定義

関数 $f(x)$ の $x$ における導関数 $f'(x)$ は以下のように定義されます：

$f'(x)=\lim _{h\to 0}{\frac {f(x+h)-f(x)}{h}}$

偏微分

多変数関数 $f(x_{1},x_{2},\dots ,x_{n})$ に対する $x_{i}$ についての偏微分は、他の変数を定数とみなして行う微分です：

${\frac {\partial f}{\partial x_{i}}}=\lim _{h\to 0}{\frac {f(x_{1},\dots ,x_{i}+h,\dots ,x_{n})-f(x_{1},\dots ,x_{i},\dots ,x_{n})}{h}}$

実践例：

using ForwardDiff

# 数値微分の例
f(x) = x^2 * sin(x)

# ForwardDiffを使用した数値微分
x = 2.0
dx = ForwardDiff.derivative(f, x)
println("f'($x) ≈ $dx")

# 中心差分法による独自実装
numerical_derivative(f, x, h=1e-6) = (f(x + h) - f(x - h)) / (2h)

dx_custom = numerical_derivative(f, x)
println("独自実装による f'($x) ≈ $dx_custom")

1.2.2 積分の基礎

積分は、関数の累積的な効果や面積、体積を計算する際に用いられる基本的な演算です。

定積分の定義

関数 $f(x)$ の区間 $[a,b]$ における定積分は以下のように定義されます：

$\int _{a}^{b}f(x)\,dx=\lim _{n\to \infty }\sum _{i=1}^{n}f(x_{i})\Delta x$

ここで、 $\Delta x={\frac {b-a}{n}}$ であり、 $x_{i}$ は各小区間の代表点です。

数値積分の実装

function trapezoidal_integration(f, a, b, n=1000)
    """台形則による数値積分"""
    x = range(a, stop=b, length=n)
    dx = (b - a) / (n - 1)
    return dx * (sum(f.(x)) - (f(a) + f(b))/2)
end

# 例：sin(x)の0からπまでの積分（結果は2になるはず）
g(x) = sin(x)

result = trapezoidal_integration(g, 0, π)
println("∫[0→π] sin(x)dx ≈ $result")

# より高度な方法：QuadGKの使用
using QuadGK
result_quadgk, error = quadgk(g, 0, π)
println("QuadGK使用: ∫[0→π] sin(x)dx ≈ $result_quadgk")

1.2.3 応用例：最適化問題

微分積分学の重要な応用例として、最適化問題を考えましょう。特に、勾配降下法は機械学習でも頻繁に使用される手法です。

function 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 1:max_iter
        gradient = df(x)
        x_new = x - learning_rate * gradient

        if abs(f(x_new) - f(x)) < tol
            break
        end

        x = x_new
        push!(history, x)
    end

    return x, history
end

# 例：x²の最小値を求める
f(x) = x^2
df(x) = 2x

x_min, history = gradient_descent(f, df, 2.0)
println("x²の最小値は x = $x_min で達成されます")

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₁)

using LinearAlgebra

# ベクトル演算の例
a = [1, 2, 3]
b = [4, 5, 6]

# 内積
dot_product = dot(a, b)
println("内積: $dot_product")

# ノルム（L2ノルム）
norm_a = norm(a)
println("ベクトルaのノルム: $norm_a")

# 外積
cross_product = cross(a, b)
println("外積: $cross_product")

行列の基本演算

行列の和差：(A±B)ᵢⱼ = Aᵢⱼ±Bᵢⱼ
行列の積：(AB)ᵢⱼ = Σₖ AᵢₖBₖⱼ
転置：(Aᵀ)ᵢⱼ = Aⱼᵢ

# 行列演算の例
A = [1 2; 3 4]
B = [5 6; 7 8]

# 行列の和
C = A + B
println("行列の和:\n", C)

# 行列の積
D = A * B
println("行列の積:\n", D)

# 転置
AT = transpose(A)
println("Aの転置:\n", AT)

1.3.2 線形方程式系の解法

線形方程式系 $Ax=b$ の解法は、科学技術計算における基本的な問題の一つです。

ガウス消去法

function gaussian_elimination(A, b)
    """
    ガウス消去法による線形方程式系の解法
    A: n×n係数行列
    b: n次元ベクトル
    """
    n = size(A, 1)
    Ab = [A b]  # 拡大係数行列

    # 前進消去
    for i in 1:n
        # ピボット選択
        pivot = argmax(abs.(Ab[i:n, i])) + i - 1
        Ab[i, :], Ab[pivot, :] = Ab[pivot, :], Ab[i, :]
        for j in i+1:n
            factor = Ab[j, i] / Ab[i, i]
            Ab[j, :] -= factor * Ab[i, :]
        end
    end

    # 後退代入
    x = zeros(n)
    for i in n:-1:1
        x[i] = (Ab[i, end] - sum(Ab[i, i+1:n] .* x[i+1:n])) / Ab[i, i]
    end

    return x
end

# 例題
A = [2 1 -1; -3 -1 2; -2 1 2]
b = [8, -11, -3]
x = gaussian_elimination(A, b)
println("解: ", x)

1.3.3 固有値と固有ベクトル

固有値問題は、振動解析や主成分分析など、多くの応用場面で現れます。

function power_method(A; max_iter=1000, tol=1e-10)
    """
    べき乗法による最大固有値と対応する固有ベクトルの計算
    """
    n = size(A, 1)
    x = rand(n)
    x /= norm(x)
    
    for _ in 1:max_iter
        x_new = A * x
        x_new /= norm(x_new)
        if norm(x - x_new) < tol
            break
        end
        x = x_new
    end

    λ = (x' * A * x) / (x' * x)
    return λ, x
end

# 例題
A = [4 -2; -2 5]
eval_power, evec_power = power_method(A)
println("べき乗法による最大固有値: ", eval_power)
println("対応する固有ベクトル: ", evec_power)

1.3.4 応用例：画像の主成分分析

線形代数の重要な応用例として、画像の次元削減を行う主成分分析を実装してみましょう。

function pca_compression(image, n_components)
    """
    PCAによる画像圧縮
    image: グレースケール画像（2次元配列）
    n_components: 使用する主成分の数
    """
    mean_img = mean(image, dims=2)
    centered = image .- mean_img
    
    # 共分散行列の計算
    cov_matrix = centered * centered' / size(centered, 2)
    
    # 固有値分解
    eigenvalues, eigenvectors = eigen(cov_matrix)
    idx = sortperm(eigenvalues, rev=true)
    eigenvectors = eigenvectors[:, idx]
    
    # 次元削減
    components = eigenvectors[:, 1:n_components]
    projected = components' * centered
    reconstructed = components * projected .+ mean_img
    
    return reconstructed
end

# 例：ランダムなパターンの圧縮
using Random
Random.seed!(42)
original = rand(100, 1000)  # 100×1000のランダム画像
compressed = pca_compression(original, 10)

# 圧縮率と再構成誤差の計算
compression_ratio = length(original) / (length(compressed) + size(compressed, 2))
error = mean((original - compressed) .^ 2)
println("圧縮率: ", compression_ratio)
println("平均二乗誤差: ", error)

1.3.2 線形方程式系の解法

線形方程式系 $Ax=b$ の解法は、科学技術計算における基本的な問題の一つです。

ガウス消去法

function gaussian_elimination(A, b)
    """
    ガウス消去法による線形方程式系の解法
    A: n×n係数行列
    b: n次元ベクトル
    """
    n = size(A, 1)
    # 拡大係数行列を作成
    Ab = hcat(A, b)

    # 前進消去
    for i in 1:n
        # ピボット選択
        pivot = argmax(abs.(Ab[i:end, i])) + i - 1
        if pivot != i
            Ab[[i, pivot], :] = Ab[[pivot, i], :]
        end

        for j in (i + 1):n
            factor = Ab[j, i] / Ab[i, i]
            Ab[j, :] -= factor * Ab[i, :]
        end
    end

    # 後退代入
    x = zeros(n)
    for i in n:-1:1
        x[i] = (Ab[i, end] - dot(Ab[i, i+1:end-1], x[i+1:end])) / Ab[i, i]
    end

    return x
end

# 例題
A = [2.0 1.0 -1.0;
     -3.0 -1.0 2.0;
     -2.0 1.0 2.0]
b = [8.0, -11.0, -3.0]

x = gaussian_elimination(A, b)
println("解: ", x)

# 検証：Juliaのバックスラッシュ演算子を使用
x_julia = A \ b
println("Juliaのバックスラッシュ演算子による解: ", x_julia)

1.3.3 固有値と固有ベクトル

固有値問題は、振動解析や主成分分析など、多くの応用場面で現れます。

べき乗法

function power_method(A, max_iter=1000, tol=1e-10)
    """
    べき乗法による最大固有値と対応する固有ベクトルの計算
    """
    n = size(A, 1)
    x = rand(n)
    x = x / norm(x)

    for _ in 1:max_iter
        x_new = A * x
        x_new = x_new / norm(x_new)

        if all(isapprox.(x, x_new, rtol=tol))
            break
        end
        x = x_new
    end

    eigenvalue = dot(x, A * x) / dot(x, x)
    return eigenvalue, x
end

# 例：対称行列の固有値・固有ベクトル
A = [4.0 -2.0;
     -2.0 5.0]

# べき乗法による計算
eval_power, evec_power = power_method(A)
println("べき乗法による最大固有値: ", eval_power)
println("対応する固有ベクトル: ", evec_power)

# Juliaの固有値分解による計算（検証用）
evals, evecs = eigen(A)
println("\nJuliaのeigen関数による固有値: ", evals)
println("対応する固有ベクトル: ", evecs)

1.3.4 応用例：画像の主成分分析

線形代数の重要な応用例として、画像の次元削減を行う主成分分析を実装してみましょう。

function pca_compression(image, n_components)
    """
    PCAによる画像圧縮
    image: グレースケール画像（2次元配列）
    n_components: 使用する主成分の数
    """
    # 平均を引く
    mean = mean(image, dims=2)
    centered = image .- mean

    # 共分散行列の計算
    cov = centered * centered' / size(centered, 2)

    # 固有値分解
    eigenvalues, eigenvectors = eigen(cov)
    idx = sortperm(eigenvalues, rev=true)
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]

    # 次元削減
    components = eigenvectors[:, 1:n_components]
    projected = components' * centered
    reconstructed = components * projected .+ mean

    return reconstructed
end

# 例：ランダムなパターンの圧縮
using Random
Random.seed!(42)
original = rand(100, 1000)  # 100×1000のランダム画像
compressed = pca_compression(original, 10)

# 圧縮率と再構成誤差の計算
compression_ratio = length(original) / (length(compressed) + size(compressed, 2))
error = mean((original .- compressed).^2)
println("圧縮率: ", compression_ratio)
println("平均二乗誤差: ", error)

1.4 複素解析

複素解析は、科学技術計算において信号処理、制御理論、量子力学など、様々な分野で重要な役割を果たします。

1.4.1 複素数の基礎

# 複素数の生成と基本演算
z1 = 2 + 3im
z2 = 1 - 2im

# 四則演算
println("和: ", z1 + z2)
println("差: ", z1 - z2)
println("積: ", z1 * z2)
println("商: ", z1 / z2)

# 極形式への変換
r = abs(z1)      # 絶対値（モジュラス）
theta = angle(z1)  # 偏角（アーギュメント）
println("絶対値: ", r)
println("偏角: ", theta, " rad")

# 極形式から直交形式への変換
z3 = r * (cos(theta) + im * sin(theta))
println("極形式から再構成: ", z3)

1.4.2 複素関数と解析性

function 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)

    println("∂u/∂x ≈ ", du_dx, ", ∂v/∂y ≈ ", dv_dy)
    println("∂u/∂y ≈ ", du_dy, ", -∂v/∂x ≈ ", -dv_dx)

    return isapprox(du_dx, dv_dy) && isapprox(du_dy, -dv_dx)
end

# 例：f(z) = z² = (x + yi)² の場合
u(x, y) = x^2 - y^2  # 実部
v(x, y) = 2*x*y      # 虚部

# z = 1 + i での解析性を確認
is_analytic = check_analytic(u, v, 1.0, 1.0)
println("関数は解析的か: ", is_analytic)

1.4.3 複素積分

function residue_theorem_example()
    """
    留数定理を用いた積分の例：
    ∮ (1/(z²+1)) dz を単位円周上で計算
    """
    # 極点: z = ±i
    # 単位円内の極: z = i
    # 留数: 1/(2i)

    using QuadGK

    integrand(theta) = begin
        # z = e^(iθ) とパラメータ化
        z = exp(im * theta)
        return 1 / (z^2 + 1) * im * z
    end

    # 数値積分（0 から 2π）
    result_numerical, _ = quadgk(integrand, 0, 2π)

    # 留数定理による結果
    result_residue = 2π * im * (1/(2im))

    println("数値積分結果: ", real(result_numerical))
    println("留数定理による結果: ", real(result_residue))
end

residue_theorem_example()

1.4.4 応用例：FFTを用いた信号処理

function signal_processing_example()
    """
    FFTを用いた信号のフィルタリング例
    """
    # サンプリングパラメータ
    t = LinRange(0, 1, 1000)
    dt = t[2] - t[1]
    freq = fftfreq(length(t), 1/dt)

    # 元の信号（2つの周波数成分）
    signal = sin.(2π*10*t) + 0.5*sin.(2π*50*t)

    # ノイズ追加
    noisy_signal = signal + randn(length(t)) * 0.2

    # FFTの実行
    fft_result = fft(noisy_signal)

    # ローパスフィルタの適用
    cutoff_freq = 30
    mask = abs.(freq) .< cutoff_freq
    filtered_fft = fft_result .* mask

    # 逆FFTで信号を復元
    filtered_signal = real(ifft(filtered_fft))

    return t, signal, noisy_signal, filtered_signal
end

# 信号処理の実行と結果の表示
t, orig, noisy, filtered = signal_processing_example()
println("信号処理完了")

1.5 確率・統計

確率・統計は、データ分析、機械学習、実験データの解析など、幅広い応用場面で必要となる基礎知識です。

1.5.1 確率の基礎

using Distributions

# サイコロを振る実験のシミュレーション
function dice_experiment(n_trials=10000)
    """
    公平なサイコロを振る実験のシミュレーション
    """
    results = rand(1:6, n_trials)
    probabilities = counts(results, 1:6) / n_trials

    println("実験的確率:")
    for face in 1:6
        println("目$face: ", probabilities[face])
    end

    # カイ二乗検定で一様性を確認
    chi2_stat = sum((counts(results, 1:6) .- n_trials/6).^2 / (n_trials/6))
    p_value = ccdf(Chisq(5), chi2_stat)

    println("\nカイ二乗検定:")
    println("統計量: ", chi2_stat)
    println("p値: ", p_value)
end

dice_experiment()

1.5.2 確率分布と統計量

function distribution_examples()
    """
    代表的な確率分布のサンプリングと可視化
    """
    using Random
    Random.seed!(42)
    n_samples = 1000

    # 正規分布
    normal_samples = rand(Normal(0, 1), n_samples)

    # 指数分布
    exponential_samples = rand(Exponential(1), n_samples)

    # 二項分布
    binomial_samples = rand(Binomial(10, 0.3), n_samples)

    # 基本統計量の計算
    function compute_statistics(samples, name)
        println("\n$nameの統計量:")
        println("平均: ", mean(samples))
        println("分散: ", var(samples))
        println("標準偏差: ", std(samples))
        println("中央値: ", median(samples))
    end

    compute_statistics(normal_samples, "正規分布")
    compute_statistics(exponential_samples, "指数分布")
    compute_statistics(binomial_samples, "二項分布")
end

distribution_examples()

1.5.3 統計的推定と検定

function statistical_inference_example()
    """
    統計的推定と検定の例
    """
    using Distributions
    using HypothesisTests

    # データ生成
    Random.seed!(42)
    population_mean = 100
    population_std = 15
    sample_size = 30

    sample = rand(Normal(population_mean, population_std), sample_size)

    # 平均の信頼区間
    confidence_level = 0.95
    t_dist = TDist(sample_size - 1)
    t_stat = quantile(t_dist, (1 + confidence_level) / 2)
    sample_mean = mean(sample)
    sample_std = std(sample)
    margin_error = t_stat * (sample_std / sqrt(sample_size))

    ci_lower = sample_mean - margin_error
    ci_upper = sample_mean + margin_error

    println("$confidence_level%信頼区間:")
    println("[", ci_lower, ", ", ci_upper, "]")

    # 仮説検定
    null_mean = 95  # 帰無仮説の平均値
    test_result = OneSampleTTest(sample, null_mean)

    println("\n一標本のt検定:")
    println("帰無仮説: μ = ", null_mean)
    println("t統計量: ", test_result.t)
    println("p値: ", pvalue(test_result))
end

statistical_inference_example()

1.5.4 回帰分析

using GLM
using DataFrames

function regression_analysis_example()
    """
    回帰分析の実装例
    """
    # データ生成
    Random.seed!(42)
    n_samples = 100

    # 説明変数
    X1 = rand(n_samples)
    X2 = rand(n_samples)

    # 目的変数（ノイズを含む）
    y = 2*X1 + 3*X2 + randn(n_samples) * 0.1

    # データフレームの作成
    data = DataFrame(X1=X1, X2=X2, y=y)

    # モデルの学習
    model = lm(@formula(y ~ X1 + X2), data)

    # 結果の表示
    println("回帰係数:")
    println(coef(model))

    println("\n決定係数 R²: ", r2(model))
    println("平均二乗誤差: ", deviance(model) / n_samples)
end

regression_analysis_example()

1.5.5 応用例：ブートストラップ法による信頼区間推定

using Random
using Statistics

function bootstrap_example()
    """
    ブートストラップ法による中央値の信頼区間推定
    """
    # データ生成
    Random.seed!(42)
    data = randn(100)
    data = exp.(data) # 対数正規分布に従う乱数

    # ブートストラップサンプリング
    n_bootstrap = 1000
    bootstrap_medians = zeros(n_bootstrap)

    for i in 1:n_bootstrap
        bootstrap_sample = rand(data, length(data)) # 重複ありランダムサンプリング
        bootstrap_medians[i] = median(bootstrap_sample)
    end

    # 信頼区間の計算（95%）
    ci_lower = quantile(bootstrap_medians, 0.025)
    ci_upper = quantile(bootstrap_medians, 0.975)

    println("中央値のブートストラップ信頼区間:")
    println("点推定値: ", median(data))
    println("95%信頼区間: [", ci_lower, ", ", ci_upper, "]")
end

bootstrap_example()

第2章数値計算の基礎

2.1 数値表現と誤差

計算機における数値計算では、有限の精度で実数を表現する必要があり、そこから生じる誤差を適切に理解し管理することが重要です。

2.1.1 浮動小数点数の表現

IEEE 754規格による浮動小数点数の表現について理解しましょう。

using Printf

function floating_point_example()
    """
    浮動小数点数の特性を確認する例
    """
    println("機械イプシロン: ", eps(Float64))
    println("最小の正規化数: ", typemin(Float64))
    println("最大の有限数: ", typemax(Float64))

    # 丸め誤差の例
    x = 0.1 + 0.2
    println("\n0.1 + 0.2 = ", x)
    println("0.3との差: ", x - 0.3)

    # 非正規化数の例
    x = Float64(1.0)
    for _ in 1:1074
        x /= 2
        @printf("x = %.2e\n", x)
    end
    @printf("最小の非正規化数: %.2e\n", x)
end

floating_point_example()

2.1.2 数値誤差の種類と対策

function error_analysis_example()
    """
    数値計算における主な誤差の例と対策
    """
    # 桁落ち（相対的に大きな数からの小さな数の減算）
    function catastrophic_cancellation()
        x = 1.0
        delta = 1e-16
        bad_result = (x + delta) - x
        println("桁落ちの例:")
        println("期待値: $delta")
        println("実際の結果: $bad_result")

        # 対策: 数式の変形
        better_result = delta
        println("改善後の結果: $better_result")
    end

    # 丸め誤差の蓄積
    function error_accumulation()
        sum1 = sum([0.1 for _ in 1:10])
        sum2 = 0.1 * 10
        println("\n丸め誤差の蓄積:")
        println("繰り返し加算: $sum1")
        println("乗算: $sum2")
        println("差: $(sum1 - sum2)")

        # 対策: Kahan summation algorithm
        function 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
            end
            return s
        end

        sum3 = kahan_sum([0.1 for _ in 1:10])
        println("Kahan summation結果: $sum3")
    end

    catastrophic_cancellation()
    error_accumulation()
end

error_analysis_example()

2.1.3 条件数と数値安定性

function condition_number_example()
    """
    条件数と数値安定性の関係を示す例
    """
    # 行列の条件数と連立方程式の解の精度
    function matrix_condition()
        # 良条件の行列
        A1 = [2.0 1.0; 1.0 2.0]

        # 悪条件の行列
        A2 = [2.0 1.999; 1.999 2.0]

        b = [1.0, 1.0]

        println("条件数の比較:")
        println("良条件行列の条件数: $(cond(A1))")
        println("悪条件行列の条件数: $(cond(A2))")

        # 摂動を加えた場合の解の変化
        delta = 1e-10
        b_perturbed = b .+ delta

        x1 = A1 \ b
        x1_perturbed = A1 \ b_perturbed

        x2 = A2 \ b
        x2_perturbed = A2 \ b_perturbed

        println("\n摂動に対する解の変化:")
        println("良条件行列:")
        println("相対誤差: $(norm(x1_perturbed - x1) / norm(x1))")
        println("悪条件行列:")
        println("相対誤差: $(norm(x2_perturbed - x2) / norm(x2))")
    end

    matrix_condition()
end

condition_number_example()

2.1.4 実践的な対策と推奨事項

function numerical_best_practices()
    """
    数値計算における推奨プラクティス
    """
    # 1. スケーリング
    function scaling_example()
        # 悪い例
        x = [1e-8, 1e8]
        println("スケーリング前の数値範囲: $(minimum(x)) to $(maximum(x))")

        # 良い例
        x_scaled = x ./ maximum(abs.(x))
        println("スケーリング後の数値範囲: $(minimum(x_scaled)) to $(maximum(x_scaled))")
    end

    # 2. 数値的に安定なアルゴリズムの選択
    function stable_algorithm_example()
        # 悪い例：不安定な二次方程式の解法
        function unstable_quadratic(a, b, c)
            x1 = (-b + sqrt(b^2 - 4*a*c)) / (2*a)
            x2 = (-b - sqrt(b^2 - 4*a*c)) / (2*a)
            return x1, x2
        end

        # 良い例：数値的に安定な解法
        function stable_quadratic(a, b, c)
            if b >= 0
                x1 = (-b - sqrt(b^2 - 4*a*c)) / (2*a)
                x2 = c / (a*x1)
            else
                x1 = (-b + sqrt(b^2 - 4*a*c)) / (2*a)
                x2 = c / (a*x1)
            end
            return x1, x2
        end

        # 比較
        a, b, c = 1.0, 1e10, 1.0
        println("\n二次方程式の解法比較:")
        println("不安定な方法: $(unstable_quadratic(a, b, c))")
        println("安定な方法: $(stable_quadratic(a, b, c))")
    end

    scaling_example()
    stable_algorithm_example()
end

numerical_best_practices()

2.2 方程式の解法

2.2.1 非線形方程式の数値解法

非線形方程式 f(x) = 0 の解を求めるための代表的な手法を実装します。

function newton_method(f, df, x0; tol=1e-6, max_iter=100)
    """
    ニュートン法による非線形方程式の求解

    Parameters:
        f: 対象となる関数
        df: fの導関数
        x0: 初期値
        tol: 収束判定の閾値
        max_iter: 最大反復回数

    Returns:
        (解, 反復回数)
    """
    x = x0
    for i in 1:max_iter
        fx = f(x)
        if abs(fx) < tol
            return x, i
        end

        dfx = df(x)
        if abs(dfx) < 1e-10  # 導関数がゼロに近い場合
            error("導関数がゼロに近いため計算を中断します")
        end

        x_new = x - fx/dfx
        if abs(x_new - x) < tol
            return x_new, i + 1
        end
        x = x_new
    end

    error("最大反復回数に達しました")
end

function secant_method(f, x0, x1; tol=1e-6, max_iter=100)
    """
    割線法による非線形方程式の求解
    """
    for i in 1:max_iter
        f0, f1 = f(x0), f(x1)
        if abs(f1) < tol
            return x1, i
        end

        if abs(f1 - f0) < 1e-10
            error("割線の傾きがゼロに近いため計算を中断します")
        end

        x_new = x1 - f1 * (x1 - x0)/(f1 - f0)
        if abs(x_new - x1) < tol
            return x_new, i + 1
        end

        x0, x1 = x1, x_new
    end

    error("最大反復回数に達しました")
end

# 方程式を解く例
function equation_solving_example()
    """
    非線形方程式を解く例
    """
    # 例1: x³ - x - 1 = 0
    f1(x) = x^3 - x - 1
    df1(x) = 3*x^2 - 1

    try
        solution_newton, iterations = newton_method(f1, df1, 1.0)
        println("ニュートン法による解:")
        println("x ≈ $solution_newton")
        println("反復回数: $iterations")
        println("残差: $(abs(f1(solution_newton)))")
    catch e
        println("エラー: $e")
    end

    # 例2: cos(x) = x
    f2(x) = cos(x) - x

    try
        solution_secant, iterations = secant_method(f2, 0.0, 1.0)
        println("\n割線法による解:")
        println("x ≈ $solution_secant")
        println("反復回数: $iterations")
        println("残差: $(abs(f2(solution_secant)))")
    catch e
        println("エラー: $e")
    end
end

equation_solving_example()

2.2.2 連立方程式の解法

function linear_system_solvers()
    """
    連立一次方程式の各種解法の実装と比較
    """
    function gauss_seidel(A, b, x0; tol=1e-6, max_iter=1000)
        """
        ガウス・ザイデル法による連立一次方程式の求解
        """
        n = length(b)
        x = copy(x0)

        for it in 1:max_iter
            x_new = copy(x)
            for i in 1:n
                s1 = dot(A[i, 1:i-1], x_new[1:i-1])
                s2 = dot(A[i, i+1:end], x[i+1:end])
                x_new[i] = (b[i] - s1 - s2) / A[i, i]
            end

            if all(abs.(x - x_new) .< tol)
                return x_new, it + 1
            end
            x = x_new
        end

        error("最大反復回数に達しました")
    end

    # 例題: 3×3の連立方程式
    A = [4.0 -1.0 0.0; -1.0 4.0 -1.0; 0.0 -1.0 4.0]
    b = [1.0, 5.0, 0.0]
    x0 = zeros(length(b))

    # 1. ガウス・ザイデル法
    try
        solution_gs, iterations = gauss_seidel(A, b, x0)
        println("ガウス・ザイデル法による解:")
        println("x = $solution_gs")
        println("反復回数: $iterations")
        println("残差ノルム: $(norm(A * solution_gs - b))")
    catch e
        println("エラー: $e")
    end

    # 2. Juliaの直接解法との比較
    solution_julia = A \ b
    println("\nJulia線形ソルバーによる解:")
    println("x = $solution_julia")
    println("残差ノルム: $(norm(A * solution_julia - b))")

    # 3. 条件数と精度の関係
    cond_A = cond(A)
    println("\n行列の条件数: $cond_A")
    println("予想される相対誤差の上限: $(cond_A * eps(Float64))")
end

linear_system_solvers()

2.2.3 最適化問題の解法

function optimization_examples()
    """
    最適化問題の解法例
    """
    using Optim

    # 1. 最急降下法の実装
    function gradient_descent(f, grad_f, x0; learning_rate=0.1, tol=1e-6, max_iter=1000)
        """
        最急降下法による最適化
        """
        x = copy(x0)
        f_value = f(x)

        for i in 1:max_iter
            grad = grad_f(x)
            if norm(grad) < tol
                return x, f_value
            end

            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
            end

            x = x_new
            f_value = f_value_new
        end

        error("最大反復回数に達しました")
    end

    # 例題1: ロゼンブロック関数の最小化
    function rosenbrock(x)
        return (1 - x[1])^2 + 100 * (x[2] - x[1]^2)^2
    end

    function rosenbrock_grad(x)
        return [-2*(1 - x[1]) - 400*x[1]*(x[2] - x[1]^2),
                200*(x[2] - x[1]^2)]
    end

    x0 = [-1.0, 1.0]

    # 最急降下法による解法
    try
        solution_gd, min_value = gradient_descent(rosenbrock, rosenbrock_grad, x0)
        println("最急降下法による最適化結果:")
        println("x = $solution_gd")
        println("f(x) = $min_value")
    catch e
        println("エラー: $e")
    end

    # Optim.jlの最適化ソルバーとの比較
    result = optimize(rosenbrock, x0, BFGS(), autodiff=:forward)
    println("\nOptim.jl BFGSによる最適化結果:")
    println("x = $(result.minimizer)")
    println("f(x) = $(result.minimum)")
    println("収束: $(Optim.converged(result))")
    println("反復回数: $(Optim.iterations(result))")
end

optimization_examples()

2.3 数値微分・積分

2.3.1 数値微分の手法

struct DerivativeResult
    value::Float64
    error_estimate::Float64
    step_size::Float64
end

function numerical_derivative(f::Function, x::Float64; method::String="central")
    eps = eps(Float64)
    h = sqrt(eps) * (1 + abs(x))  # 初期ステップサイズ
    
    if method == "forward"
        derivative = (f(x + h) - f(x)) / h
        error_est = abs(h)
        
    elseif method == "central"
        derivative = (f(x + h) - f(x - h)) / (2 * h)
        error_est = abs(h)^2
        
    elseif 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
        error("Unknown method")
    end
    
    return DerivativeResult(derivative, error_est, h)
end

function higher_derivatives(f::Function, x::Float64, order::Int=2)
    if order < 1
        error("Order must be positive")
    end
    
    h = cbrt(eps(Float64)) * (1 + abs(x))
    
    if order == 1
        return numerical_derivative(f, x, method="central").value
    elseif order == 2
        # 中心差分による2階導関数
        return (f(x + h) - 2*f(x) + f(x - h)) / (h * h)
    else
        # 再帰的に高次導関数を計算
        first_derivative(y) = numerical_derivative(f, y, method="central").value
        return higher_derivatives(first_derivative, x, order - 1)
    end
end

# 数値微分の例
function derivative_examples()
    # テスト関数
    f(x) = sin(x) * exp(-x/2)
    # 解析的な導関数
    df(x) = exp(-x/2) * (cos(x) - 0.5*sin(x))
    
    x = 1.0
    
    # 各手法での計算結果
    methods = ["forward", "central", "richardson"]
    println("数値微分の比較:")
    for method in methods
        result = numerical_derivative(f, x, method=method)
        println("\n$(method)差分法:")
        println("計算値: $(result.value)")
        println("真値との誤差: $(abs(result.value - df(x)))")
        println("推定誤差: $(result.error_estimate)")
        println("使用ステップサイズ: $(result.step_size)")
    end
    
    # 高次導関数の計算
    println("\n高次導関数:")
    for order in 1:3
        deriv = higher_derivatives(f, x, order)
        println("$(order)階導関数: $(deriv)")
    end
end

derivative_examples()

2.3.2 数値積分の手法

struct IntegrationResult
    value::Float64
    error_estimate::Float64
    n_evaluations::Int
end

function adaptive_simpson(f::Function, a::Float64, b::Float64; tol::Float64=1e-8, max_depth::Int=50)
    function simpson_rule(f, a, b)
        c = (a + b) / 2
        h = (b - a) / 6
        return h * (f(a) + 4*f(c) + f(b))
    end
    
    function _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)
        end
        
        if abs(whole_integral - integral) <= 15 * tol
            return whole_integral, abs(whole_integral - integral) / 15
        end
        
        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
    end
    
    # 初期評価
    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)
end

function gauss_quadrature(f::Function, a::Float64, b::Float64; n::Int=5)
    # ガウス・ルジャンドル求積の重みと節点（5点法）
    x5 = [0.0, ±0.5384693101056831, ±0.9061798459386640]
    w5 = [0.5688888888888889, 0.4786286704993665, 0.2369268850561891]
    
    # 区間変換
    mid = (b + a) / 2
    half_length = (b - a) / 2
    x = mid .+ half_length .* x5
    
    # 積分計算
    result = half_length * sum(w5 .* f.(x))
    
    # 誤差推定（簡易的）
    error_est = abs(result) * 1e-12
    
    return IntegrationResult(result, error_est, 5)
end

# 数値積分の例
function integration_examples()
    # テスト関数
    f(x) = exp(-x^2) * sin(2x)
    
    # 積分区間
    a, b = 0, 2
    
    # 適応的シンプソン法
    result_simpson = adaptive_simpson(f, a, b)
    println("適応的シンプソン法:")
    println("積分値: $(result_simpson.value)")
    println("推定誤差: $(result_simpson.error_estimate)")
    
    # ガウス求積法
    result_gauss = gauss_quadrature(f, a, b)
    println("\nガウス求積法:")
    println("積分値: $(result_gauss.value)")
    println("推定誤差: $(result_gauss.error_estimate)")
    
    # QuadGKとの比較
    using QuadGK
    result_quadgk, error_quadgk = quadgk(f, a, b)
    println("\nQuadGK:")
    println("積分値: $(result_quadgk)")
    println("推定誤差: $(error_quadgk)")
end

integration_examples()

2.4 線形方程式の解法

2.4.1 直接法

function lu_decomposition(A::Matrix{Float64})
    n = size(A, 1)
    L = zeros(n, n)
    U = zeros(n, n)
    
    for i in 1:n
        # L の対角成分
        L[i, i] = 1.0
        
        # U の i 行目
        for j in i:n
            sum_u = sum(L[i, k] * U[k, j] for k in 1:i-1)
            U[i, j] = A[i, j] - sum_u
        end
        
        # L の i 列目
        for j in i+1:n
            sum_l = sum(L[j, k] * U[k, i] for k in 1:i-1)
            L[j, i] = (A[j, i] - sum_l) / U[i, i]
        end
    end
    
    return L, U
end

function solve_lu(L::Matrix{Float64}, U::Matrix{Float64}, b::Vector{Float64})
    n = length(b)
    # 前進代入 (Ly = b)
    y = zeros(n)
    for i in 1:n
        y[i] = b[i] - dot(L[i, 1:i-1], y[1:i-1])
    end
    
    # 後退代入 (Ux = y)
    x = zeros(n)
    for i in n:-1:1
        x[i] = (y[i] - dot(U[i, i+1:n], x[i+1:n])) / U[i, i]
    end
    
    return x
end

2.4.2 反復法

function conjugate_gradient(A::Matrix{Float64}, b::Vector{Float64}; x0::Union{Vector{Float64}, Nothing}=nothing, tol::Float64=1e-10, max_iter::Int=1000)
    n = length(b)
    if x0 === nothing
        x = zeros(n)
    else
        x = copy(x0)
    end
    
    r = b - A * x
    p = copy(r)
    r_norm_sq = dot(r, r)
    
    for i in 1:max_iter
        Ap = A * p
        alpha = r_norm_sq / dot(p, Ap)
        x += alpha * p
        r -= alpha * Ap
        
        r_norm_sq_new = dot(r, r)
        if sqrt(r_norm_sq_new) < tol
            return x, i + 1
        end
        
        beta = r_norm_sq_new / r_norm_sq
        r_norm_sq = r_norm_sq_new
        p = r + beta * p
    end
    
    error("最大反復回数に達しました")
end

2.4.3 前処理付き反復法

function incomplete_cholesky(A::Matrix{Float64})
    n = size(A, 1)
    L = zeros(n, n)
    
    for i in 1:n
        for j in 1:i
            s = A[i, j]
            for k in 1:j-1
                s -= L[i, k] * L[j, k]
            end
            if i == j
                if s <= 0
                    error("行列が正定値でない可能性があります")
                end
                L[i, j] = sqrt(s)
            else
                L[i, j] = s / L[j, j]
            end
        end
    end
    
    return L
end

function preconditioned_cg(A::Matrix{Float64}, b::Vector{Float64}, M::Matrix{Float64}; x0::Union{Vector{Float64}, Nothing}=nothing, tol::Float64=1e-10, max_iter::Int=1000)
    n = length(b)
    if x0 === nothing
        x = zeros(n)
    else
        x = copy(x0)
    end
    
    r = b - A * x
    z = M \ r
    p = copy(z)
    rz = dot(r, z)
    
    for i in 1:max_iter
        Ap = A * p
        alpha = rz / dot(p, Ap)
        x += alpha * p
        r -= alpha * Ap
        
        if norm(r) < tol
            return x, i + 1
        end
        
        z = M \ r
        rz_new = dot(r, z)
        beta = rz_new / rz
        rz = rz_new
        p = z + beta * p
    end
    
    error("最大反復回数に達しました")
end

2.5 固有値計算

2.5.1 べき乗法

function power_method(A::Matrix{Float64}; tol::Float64=1e-10, max_iter::Int=1000)
    n = size(A, 1)
    x = rand(n)
    x = x / norm(x)
    
    lambda_old = 0
    
    for i in 1:max_iter
        # 行列とベクトルの積
        y = A * x
        
        # レイリー商による固有値の計算
        lambda_new = dot(x, y)
        
        # ベクトルの正規化
        x = y / norm(y)
        
        # 収束判定
        if abs(lambda_new - lambda_old) < tol
            return lambda_new, x
        end
        
        lambda_old = lambda_new
    end
    
    error("最大反復回数に達しました")
end

第3章信号処理

3.1 フーリエ解析

フーリエ解析は信号処理の基礎となる重要な手法です。時間領域の信号を周波数領域で解析することで、信号の特性を理解し、効果的な処理が可能になります。

3.1.1 離散フーリエ変換（DFT）の基礎

using Plots

function dft(x::AbstractVector{ComplexF64})
    """
    離散フーリエ変換の直接計算による実装

    Parameters:
        x: 入力信号

    Returns:
        X: 周波数領域の信号
    """
    N = length(x)
    X = zeros(ComplexF64, N)

    for k in 0:N-1
        for n in 0:N-1
            X[k+1] += x[n+1] * exp(-2im * pi * k * n / N)
        end
    end

    return X
end

function idft(X::AbstractVector{ComplexF64})
    """
    逆離散フーリエ変換の直接計算による実装

    Parameters:
        X: 周波数領域の信号

    Returns:
        x: 時間領域の信号
    """
    N = length(X)
    x = zeros(ComplexF64, N)

    for n in 0:N-1
        for k in 0:N-1
            x[n+1] += X[k+1] * exp(2im * pi * k * n / N)
        end
    end

    return x / N
end

function analyze_signal(t::AbstractVector{Float64}, x::AbstractVector{Float64})
    """
    信号のスペクトル解析を行う

    Parameters:
        t: 時間配列
        x: 信号配列

    Returns:
        freq: 周波数配列
        spectrum: パワースペクトル
    """
    N = length(t)
    dt = t[2] - t[1]  # サンプリング間隔

    # FFTの計算
    X = fft(x)

    # 周波数軸の生成
    freq = fftfreq(N, 1/dt)

    # パワースペクトルの計算
    spectrum = abs2.(X)

    return freq, spectrum
end

3.1.2 窓関数と周波数漏れ

function apply_window(x::AbstractVector{Float64}, window_type::String="hanning")
    """
    信号に窓関数を適用する

    Parameters:
        x: 入力信号
        window_type: 窓関数の種類 ('hanning', 'hamming', 'blackman')

    Returns:
        x_windowed: 窓関数適用後の信号
    """
    N = length(x)

    if window_type == "hanning"
        window = 0.5 * (1 .- cos.(2 * pi * (0:N-1) / (N - 1)))
    elseif window_type == "hamming"
        window = 0.54 .- 0.46 * cos.(2 * pi * (0:N-1) / (N - 1))
    elseif window_type == "blackman"
        window = (0.42 .- 0.5 * cos.(2 * pi * (0:N-1) / (N - 1)) .+
                  0.08 * cos.(4 * pi * (0:N-1) / (N - 1)))
    else
        error("未対応の窓関数です")
    end

    return x .* window
end

function spectral_analysis_with_window(t::AbstractVector{Float64},
                                      x::AbstractVector{Float64},
                                      window_type::String="hanning")
    """
    窓関数を適用してスペクトル解析を行う

    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
end

3.1.3 短時間フーリエ変換（STFT）

function stft(x::AbstractVector{Float64},
              window_size::Int64,
              hop_size::Int64,
              window_type::String="hanning")
    """
    短時間フーリエ変換の実装

    Parameters:
        x: 入力信号
        window_size: 窓サイズ
        hop_size: フレームシフト量
        window_type: 窓関数の種類

    Returns:
        times: 時間フレーム
        frequencies: 周波数ビン
        spectrogram: スペクトログラム
    """
    # 窓関数の生成
    window = apply_window(ones(window_size), window_type)

    # フレーム数の計算
    n_frames = 1 + (length(x) - window_size) ÷ hop_size

    # スペクトログラムの初期化
    spectrogram = zeros(ComplexF64, window_size ÷ 2 + 1, n_frames)

    # STFTの計算
    for i in 0:n_frames-1
        start = i * hop_size
        frame = x[start+1:start + window_size] .* window
        spectrum = rfft(frame)
        spectrogram[:, i+1] = spectrum
    end

    # 時間軸と周波数軸の生成
    times = (0:n_frames-1) * hop_size
    frequencies = rfftfreq(window_size, 1.0)

    return times, frequencies, spectrogram
end

function istft(spectrogram::AbstractMatrix{ComplexF64},
              window_size::Int64,
              hop_size::Int64,
              window_type::String="hanning")
    """
    逆短時間フーリエ変換の実装

    Parameters:
        spectrogram: スペクトログラム
        window_size: 窓サイズ
        hop_size: フレームシフト量
        window_type: 窓関数の種類

    Returns:
        x_reconstructed: 再構成された信号
    """
    # 窓関数の生成
    window = apply_window(ones(window_size), window_type)

    n_frames = size(spectrogram, 2)
    expected_signal_len = hop_size * (n_frames - 1) + window_size

    # 信号の再構成
    x_reconstructed = zeros(expected_signal_len)
    window_sum = zeros(expected_signal_len)

    for i in 0:n_frames-1
        start = i * hop_size
        frame = irfft(spectrogram[:, i+1], window_size)
        x_reconstructed[start+1:start + window_size] += frame .* window
        window_sum[start+1:start + window_size] += window .^ 2
    end

    # オーバーラップアド処理
    non_zero = window_sum .> 1e-10
    x_reconstructed[non_zero] ./= window_sum[non_zero]

    return x_reconstructed
end

3.2 ディジタルフィルタ

3.2.1 FIRフィルタの設計

function design_fir_filter(cutoff::Float64,
                           filter_order::Int64,
                           filter_type::String="lowpass")
    """
    窓関数法によるFIRフィルタの設計

    Parameters:
        cutoff: カットオフ周波数 (0 < cutoff < 0.5)
        filter_order: フィルタの次数
        filter_type: フィルタの種類 ('lowpass' or 'highpass')

    Returns:
        h: フィルタ係数
    """
    if !(0 < cutoff < 0.5)
        error("カットオフ周波数は0と0.5の間である必要があります")
    end

    # インパルス応答の計算
    n = 0:filter_order
    omega_c = 2 * pi * cutoff

    if filter_type == "lowpass"
        h = sinc.(omega_c * ((n .- filter_order/2) / pi)
    elseif filter_type == "highpass"
        h = sinc.(n .- filter_order/2) .- sinc.(omega_c * (n .- filter_order/2) / pi)
    else
        error("未対応のフィルタ種類です")
    end

    # ハミング窓の適用
    window = 0.54 .- 0.46 * cos.(2 * pi * n / filter_order)
    h = h .* window

    # 正規化
    h = h / sum(h)

    return h
end

function apply_fir_filter(x::AbstractVector{Float64}, h::AbstractVector{Float64})
    """
    FIRフィルタの適用

    Parameters:
        x: 入力信号
        h: フィルタ係数

    Returns:
        y: フィルタリング後の信号
    """
    return conv(x, h, "same")
end

3.2.2 IIRフィルタの設計

using DSP

function design_iir_butterworth(cutoff::Float64,
                                 order::Int64,
                                 filter_type::String="lowpass")
    """
    バターワースIIRフィルタの設計

    Parameters:
        cutoff: カットオフ周波数 (0 < cutoff < 0.5)
        order: フィルタの次数
        filter_type: フィルタの種類

    Returns:
        b: 分子の係数
        a: 分母の係数
    """
    nyquist = 0.5
    cutoff_normalized = cutoff / nyquist

    if filter_type == "lowpass"
        b, a = butter(order, cutoff_normalized, lowpass=true)
    elseif filter_type == "highpass"
        b, a = butter(order, cutoff_normalized, highpass=true)
    else
        error("未対応のフィルタ種類です")
    end

    return b, a
end

function apply_iir_filter(x::AbstractVector{Float64},
                          b::AbstractVector{Float64},
                          a::AbstractVector{Float64},
                          initial_conditions::Union{AbstractVector{Float64}, Nothing}=nothing)
    """
    IIRフィルタの適用

    Parameters:
        x: 入力信号
        b: 分子の係数
        a: 分母の係数
        initial_conditions: 初期状態

    Returns:
        y: フィルタリング後の信号
    """
    if initial_conditions === nothing
        y = filt(b, a, x)
    else
        y = filt(b, a, x, zi=initial_conditions)
    end
    return y
end

3.3 画像処理の基礎

3.3.1 空間フィルタリング

using Images

function spatial_filter(image::AbstractMatrix{Float64},
                        kernel::AbstractMatrix{Float64})
    """
    2次元畳み込みによる空間フィルタリング

    Parameters:
        image: 入力画像
        kernel: フィルタカーネル

    Returns:
        filtered_image: フィルタリング後の画像
    """
    # パディングサイズの計算
    pad_h = size(kernel, 1) ÷ 2
    pad_w = size(kernel, 2) ÷ 2

    # パディング処理
    padded_image = padarray(image, pad_h, pad_w, method=:reflect)

    # 出力画像の初期化
    filtered_image = zeros(Float64, size(image))

    # 畳み込み演算
    for i in 1:size(image, 1)
        for j in 1:size(image, 2)
            window = padded_image[i:i + size(kernel, 1) - 1, j:j + size(kernel, 2) - 1]
            filtered_image[i, j] = sum(window .* kernel)
        end
    end

    return filtered_image
end

function create_gaussian_kernel(size::Int64, sigma::Float64)
    """
    ガウシアンフィルタカーネルの生成

    Parameters:
        size: カーネルサイズ
        sigma: 標準偏差

    Returns:
        kernel: ガウシアンカーネル
    """
    x, y = meshgrid(-size ÷ 2:size ÷ 2, -size ÷ 2:size ÷ 2)

    kernel = exp.(-(x.^2 + y.^2) / (2 * sigma^2))
    return kernel / sum(kernel)
end

3.3.2 エッジ検出

function sobel_edge_detection(image::AbstractMatrix{Float64})
    """
    Sobelオペレータによるエッジ検出

    Parameters:
        image: 入力画像（グレースケール）

    Returns:
        magnitude: エッジの強度
        direction: エッジの方向
        edges: 二値化されたエッジ画像
    """
    # Sobelカーネルの定義
    kernel_x = [-1 0 1;
                -2 0 2;
                -1 0 1]

    kernel_y = [-1 -2 -1;
                0 0 0;
                1 2 1]

    # 勾配の計算
    gradient_x = spatial_filter(image, kernel_x)
    gradient_y = spatial_filter(image, kernel_y)

    # エッジの強度と方向の計算
    magnitude = sqrt.(gradient_x.^2 + gradient_y.^2)
    direction = atan.(gradient_y, gradient_x)

    # 閾値処理によるエッジの二値化
    threshold = 0.5 * maximum(magnitude)
    edges = magnitude .> threshold

    return magnitude, direction, edges
end

function canny_edge_detection(image::AbstractMatrix{Float64},
                              low_threshold::Float64=0.1,
                              high_threshold::Float64=0.3)
    """
    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 * maximum(suppressed)
    low = low_threshold * maximum(suppressed)

    return hysteresis_thresholding(suppressed, low, high)
end

function non_maximum_suppression(magnitude::AbstractMatrix{Float64},
                                  direction::AbstractMatrix{Float64})
    """
    非最大値抑制の実装
    """
    height, width = size(magnitude)
    suppressed = zeros(Float64, size(magnitude))

    # 勾配方向を0°、45°、90°、135°に量子化
    direction = rad2deg.(direction) .% 180

    for i in 2:height - 1
        for j in 2:width - 1
            if direction[i, j] < 22.5 || direction[i, j] >= 157.5
                neighbors = [magnitude[i, j - 1], magnitude[i, j + 1]]
            elseif 22.5 <= direction[i, j] < 67.5
                neighbors = [magnitude[i - 1, j - 1], magnitude[i + 1, j + 1]]
            elseif 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]]
            end

            if magnitude[i, j] >= maximum(neighbors)
                suppressed[i, j] = magnitude[i, j]
            end
        end
    end

    return suppressed
end

3.3.3 特徴点検出

function harris_corner_detector(image::AbstractMatrix{Float64},
                                 k::Float64=0.04, threshold::Float64=0.01)
    """
    Harrisコーナー検出器の実装

    Parameters:
        image: 入力画像
        k: Harrisレスポンス関数のパラメータ
        threshold: コーナー検出の閾値

    Returns:
        corners: コーナー点の位置（True/False）
    """
    # 勾配の計算
    dx = spatial_filter(image, [-1 0 1])
    dy = spatial_filter(image, [-1 0 1]')

    # 構造テンソルの要素
    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 = zeros(Bool, size(image))
    response_max = maximum(response)
    threshold_abs = threshold * response_max

    for i in 2:size(image, 1) - 1
        for j in 2:size(image, 2) - 1
            if response[i, j] > threshold_abs
                window = response[i - 1:i + 1, j - 1:j + 1]
                if response[i, j] == maximum(window)
                    corners[i, j] = true
                end
            end
        end
    end

    return corners
end

3.4 時系列解析

3.4.1 自己相関分析

function autocorrelation(x::AbstractVector{Float64},
                       max_lag::Union{Int64, Nothing}=nothing)
    """
    時系列データの自己相関関数を計算

    Parameters:
        x: 入力時系列データ
        max_lag: 最大ラグ（デフォルトはデータ長の半分）

    Returns:
        acf: 自己相関関数
    """
    if max_lag === nothing
        max_lag = length(x) ÷ 2
    end

    # 平均を引く
    x_centered = x .- mean(x)

    # 自己相関の計算
    acf = zeros(Float64, max_lag + 1)
    variance = var(x_centered)

    for lag in 0:max_lag
        acf[lag + 1] = sum(x_centered[lag + 1:end] .* x_centered[1:end - lag])
    end

    # 正規化
    return acf / (length(x) * variance)
end

function partial_autocorrelation(x::AbstractVector{Float64},
                                 max_lag::Union{Int64, Nothing}=nothing)
    """
    偏自己相関関数の計算（Durbin-Levinsonアルゴリズム）

    Parameters:
        x: 入力時系列データ
        max_lag: 最大ラグ

    Returns:
        pacf: 偏自己相関関数
    """
    if max_lag === nothing
        max_lag = length(x) ÷ 2
    end

    # 自己相関関数の計算
    acf = autocorrelation(x, max_lag)

    # 偏自己相関係数の初期化
    pacf = zeros(Float64, max_lag + 1)
    pacf[1] = 1.0

    # Durbin-Levinsonアルゴリズム
    for k in 1:max_lag
        phi = zeros(Float64, k)
        phi[k] = acf[k + 1]

        for j in 1:k - 1
            phi[j] = pacf[j]
        end

        for j in 1:k - 1
            phi[j] -= pacf[k] * pacf[k - j]
        end

        pacf[k] = phi[k]
    end

    return pacf
end

3.4.2 スペクトル推定

function welch_periodogram(x::AbstractVector{Float64},
                           window_size::Union{Int64, Nothing}=nothing,
                           overlap::Float64=0.5,
                           window_type::String="hanning")
    """
    Welch法によるパワースペクトル密度の推定

    Parameters:
        x: 入力時系列データ
        window_size: 窓サイズ
        overlap: オーバーラップ率（0～1）
        window_type: 窓関数の種類

    Returns:
        frequencies: 周波数配列
        psd: パワースペクトル密度
    """
    if window_size === nothing
        window_size = length(x) ÷ 8
    end

    # オーバーラップするサンプル数
    hop_size = Int64(floor(window_size * (1 - overlap)))

    # 窓関数の適用
    window = apply_window(ones(window_size), window_type)

    # セグメントの数の計算
    n_segments = (length(x) - window_size) ÷ hop_size + 1

    # スペクトルの初期化
    spectrum = zeros(Float64, window_size ÷ 2 + 1)

    # 各セグメントのスペクトルの計算と平均化
    for i in 0:n_segments - 1
        start = i * hop_size
        segment = x[start + 1:start + window_size] .* window
        segment_spectrum = abs2.(rfft(segment))
        spectrum += segment_spectrum
    end

    # 正規化
    spectrum /= n_segments

    # 周波数軸の生成
    frequencies = rfftfreq(window_size, 1.0)

    return frequencies, spectrum
end

3.5 ウェーブレット変換

3.5.1 連続ウェーブレット変換

function morlet_wavelet(t::AbstractVector{Float64}, omega0::Float64=6.0)
    """
    モルレーウェーブレットの生成

    Parameters:
        t: 時間配列
        omega0: 中心周波数

    Returns:
        psi: ウェーブレット関数
    """
    return pi^(-0.25) * exp.(1im * omega0 * t) .* exp.(-t.^2 / 2)
end

function cwt(x::AbstractVector{Float64},
             scales::AbstractVector{Float64},
             dt::Float64=1.0)
    """
    連続ウェーブレット変換の実装

    Parameters:
        x: 入力信号
        scales: スケールパラメータの配列
        dt: サンプリング間隔

    Returns:
        W: ウェーブレット係数（時間-スケール表現）
    """
    n = length(x)
    W = zeros(ComplexF64, length(scales), n)

    # 各スケールでのウェーブレット変換
    for i in 1:length(scales)
        scale = scales[i]
        # ウェーブレットの生成
        t = (-n ÷ 2 + 1:n ÷ 2) * dt
        wavelet = morlet_wavelet(t / scale)

        # 畳み込み（時間領域で計算）
        x_pad = padarray(x, n ÷ 2, method=:reflect)
        convolution = conv(x_pad, reverse(wavelet))
        W[i, :] = convolution[n ÷ 2 + 1:end - n ÷ 2]
    end

    return W
end

下位階層のページ

はじめに

第1章 基礎数学

1.1 代数学の基礎

1.1.1 集合と写像

1.1.2 群・環・体

1.2 微分積分学

1.2.1 微分の基礎

1.2.2 積分の基礎

1.2.3 応用例：最適化問題

1.3 線形代数

1.3.1 ベクトルと行列の基礎

1.3.1 ベクトルと行列の基本演算

1.3.2 線形方程式系の解法

1.3.3 固有値と固有ベクトル

1.3.4 応用例：画像の主成分分析

1.3.2 線形方程式系の解法

ガウス消去法

1.3.3 固有値と固有ベクトル

べき乗法

1.3.4 応用例：画像の主成分分析

1.4 複素解析

1.4.1 複素数の基礎

1.4.2 複素関数と解析性

1.4.3 複素積分

1.4.4 応用例：FFTを用いた信号処理

1.5 確率・統計

1.5.1 確率の基礎

1.5.2 確率分布と統計量

1.5.3 統計的推定と検定

1.5.4 回帰分析

1.5.5 応用例：ブートストラップ法による信頼区間推定

第2章 数値計算の基礎

2.1 数値表現と誤差

2.1.1 浮動小数点数の表現

2.1.2 数値誤差の種類と対策

2.1.3 条件数と数値安定性

2.1.4 実践的な対策と推奨事項

2.2 方程式の解法

2.2.1 非線形方程式の数値解法

2.2.2 連立方程式の解法

2.2.3 最適化問題の解法

2.3 数値微分・積分

2.3.1 数値微分の手法

2.3.2 数値積分の手法

2.4 線形方程式の解法

2.4.1 直接法

2.4.2 反復法

2.4.3 前処理付き反復法

2.5 固有値計算

2.5.1 べき乗法

第3章 信号処理

3.1 フーリエ解析

3.1.1 離散フーリエ変換（DFT）の基礎

3.1.2 窓関数と周波数漏れ

3.1.3 短時間フーリエ変換（STFT）

3.2 ディジタルフィルタ

3.2.1 FIRフィルタの設計

3.2.2 IIRフィルタの設計

3.3 画像処理の基礎

3.3.1 空間フィルタリング

3.3.2 エッジ検出

3.3.3 特徴点検出

3.4 時系列解析

3.4.1 自己相関分析

3.4.2 スペクトル推定

3.5 ウェーブレット変換

3.5.1 連続ウェーブレット変換

下位階層のページ

第1章基礎数学

第2章数値計算の基礎

第3章信号処理