Rでフーリエ級数展開

library(ggplot2)
library(dplyr)
library(tidyr)
library(purrr)

square_wave <- function(x) {
  x_mod <- ((x + pi) %% (2 * pi)) - pi
  ifelse(x_mod > 0, 1, -1)
}

sawtooth_wave <- function(x) {
  x_mod <- ((x + pi) %% (2 * pi)) - pi
  x_mod / pi # -1 から 1 へ直線的に増加し、周期端で不連続に戻る
}

triangle_wave <- function(x) {
  x_mod <- ((x + pi / 2) %% (2 * pi)) - pi / 2
  ifelse(x_mod <= pi / 2,
    (2 / pi) * x_mod,
    (2 / pi) * (pi - x_mod)
  )
}

target_functions <- list(
  "矩形波" = square_wave,
  "のこぎり波" = sawtooth_wave,
  "三角波" = triangle_wave
)

bn_analytic <- function(wave_name, n) {
  switch(wave_name,
    "矩形波" = ifelse(n %% 2 == 1, 4 / (n * pi), 0),
    "のこぎり波" = 2 * (-1)^(n + 1) / (n * pi),
    "三角波" = ifelse(n %% 2 == 1,
      8 * (-1)^((n - 1) / 2) / (n^2 * pi^2),
      0
    )
  )
}

fourier_partial_sum <- function(x, N, wave_name) {
  s <- 0
  for (n in seq_len(N)) {
    b_n <- bn_analytic(wave_name, n)
    s <- s + b_n * sin(n * x)
  }
  s
}

compute_bn_numeric <- function(f, n, lower = -pi, upper = pi) {
  integrand <- function(x) f(x) * sin(n * x)
  (1 / pi) * integrate(integrand, lower, upper,
    subdivisions = 500
  )$value
}

cat("=== 各波形のフーリエ係数:解析解 vs 数値積分 ===\n")
for (wave_name in names(target_functions)) {
  cat(sprintf("\n--- %s ---\n", wave_name))
  for (n in 1:4) {
    b_a <- bn_analytic(wave_name, n)
    b_n <- compute_bn_numeric(target_functions[[wave_name]], n)
    cat(sprintf("n=%d : 解析解 = %.6f, 数値積分 = %.6f\n", n, b_a, b_n))
  }
}

# 矩形波: 不連続点は x = 0(および x = ±π)
# のこぎり波: 不連続点は x = ±π(x = 0 は連続な滑らかな点なので注意)
# 三角波: 連続関数 → オーバーシュートは N→∞ で消える

# --- 矩形波用:x = 0 の不連続点を右側から測定 ---
overshoot_check_square <- function(N) {
  x_fine <- seq(0.0001, pi / N, length.out = 300)
  y_fine <- fourier_partial_sum(x_fine, N, "矩形波")
  max(y_fine) - 1 # 真値の上限(1)からのオーバーシュート
}

# --- のこぎり波用:x = π の不連続点を左側から測定 ---
overshoot_check_sawtooth <- function(N) {
  x_fine <- seq(pi - pi / N, pi - 0.0001, length.out = 300)
  y_fine <- fourier_partial_sum(x_fine, N, "のこぎり波")
  max(y_fine) - 1 # 真値の上限(1)からのオーバーシュート
}

cat("=== 不連続点付近のオーバーシュート量 ===\n")

cat("\n--- 矩形波(不連続点: x = 0) ---\n")
for (N in c(10, 50, 200, 1000)) {
  os <- overshoot_check_square(N)
  cat(sprintf("N = %5d : オーバーシュート = %.4f\n", N, os))
}

cat("\n--- のこぎり波(不連続点: x = π) ---\n")
for (N in c(10, 50, 200, 1000)) {
  os <- overshoot_check_sawtooth(N)
  cat(sprintf("N = %5d : オーバーシュート = %.4f\n", N, os))
}

# --- 三角波:連続関数なのでオーバーシュートは生じない ---
cat("\n--- 三角波(連続関数、Gibbs現象なし) ---\n")
for (N in c(1, 3, 5, 15)) {
  x_fine <- seq(0.0001, pi / N, length.out = 300)
  y_fine <- fourier_partial_sum(x_fine, N, "三角波")
  # 三角波は x=0付近で値0から傾き2/piで単調増加、オーバーシュートはほぼ生じない
  diff_from_true <- max(y_fine - triangle_wave(x_fine))
  cat(sprintf("N = %5d : 真値との最大差分 = %.5f\n", N, diff_from_true))
}

# ジャンプの大きさが 2 の不連続点における理論上のオーバーシュート:
#   (2/pi * Si(pi) - 1) * ジャンプ幅 ≈ 0.0895 * 2 = 0.1790
gibbs_theoretical <- function(jump_size) {
  si_pi <- integrate(function(t) sin(t) / t, 0.0001, pi)$value
  (2 / pi * si_pi - 1) * (jump_size / 2) # ← ジャンプ幅を2で正規化
}

cat(sprintf(
  "理論上のオーバーシュート量(ジャンプ幅=2の場合): %.4f\n",
  gibbs_theoretical(2)
))