Rで有限要素法:ガウス求積

# ============================================================
# 5点ガウス求積の検証:解析的な積分結果との比較
# ============================================================
#
# 5点ガウス求積は、2n-1=9次以下の多項式であれば厳密に(誤差なく)積分できる
# という性質がある。これを、実際に解析的な積分結果と比較して確認する。
#
# 検証に使う関数f(ξ)として、単純すぎないよう、0次から9次までの項を
# すべて含む9次多項式を用意する。多項式ですので、各項ξ^kの積分は
# 手計算でも簡単に求まる(kが奇数なら0、偶数なら2/(k+1))ため、
# 解析解を数値的に正確に計算できる。

# ------------------------------------------------------------
# 1. 検証用の9次多項式 f(ξ) の係数を定義する
# ------------------------------------------------------------
# f(ξ) = 5 - 2ξ + 6ξ^2 - ξ^3 + 7ξ^4 - 4ξ^5 + 2ξ^6 - 5ξ^7 + 0ξ^8 + 3ξ^9
coeffs <- c(5, -2, 6, -1, 7, -4, 2, -5, 0, 3) # 0次からの係数(coeffs[k+1]がξ^kの係数)

f <- function(xi) {
  total <- 0
  for (k in 0:9) {
    total <- total + coeffs[k + 1] * xi^k
  }
  total
}

# ------------------------------------------------------------
# 2. 解析的な積分値を求める
# ------------------------------------------------------------
# ∫[-1,1] ξ^k dξ は、kが奇数なら0、偶数なら 2/(k+1) になる。
# これを各項について計算し、係数を掛けて足し合わせる。

analytic_integral <- 0
for (k in 0:9) {
  if (k %% 2 == 0) { # kが偶数の項だけが寄与する
    analytic_integral <- analytic_integral + coeffs[k + 1] * (2 / (k + 1))
  }
}
cat("解析的な積分値:", analytic_integral, "\n")

# ------------------------------------------------------------
# 3. 5点ガウス求積による近似値を求める
# ------------------------------------------------------------
gauss_pts <- c(
  0,
  -1 / 3 * sqrt(5 - 2 * sqrt(10 / 7)),
  1 / 3 * sqrt(5 - 2 * sqrt(10 / 7)),
  -1 / 3 * sqrt(5 + 2 * sqrt(10 / 7)),
  1 / 3 * sqrt(5 + 2 * sqrt(10 / 7))
)
gauss_wts <- c(
  128 / 225,
  (322 + 13 * sqrt(70)) / 900,
  (322 + 13 * sqrt(70)) / 900,
  (322 - 13 * sqrt(70)) / 900,
  (322 - 13 * sqrt(70)) / 900
)

gauss_integral <- sum(gauss_wts * f(gauss_pts)) # 評価点での関数値×重み、を全点分足し合わせる
cat("5点ガウス求積の近似値:", gauss_integral, "\n")

# ------------------------------------------------------------
# 4. 両者の差を確認する
# ------------------------------------------------------------
cat("差:", gauss_integral - analytic_integral, "\n")

# ------------------------------------------------------------
# 5. 10次の項を加えると、理論通り誤差が生じることを確認する
# ------------------------------------------------------------
# 5点ガウス求積は「9次まで厳密」ですので、10次の項(偶数次)を加えると
# 理論上、有限の誤差が生じるはずである(奇数次の項は対称性により
# 常にゼロに積分されるため、ここでは偶数次である10次の項を追加する)。

coeffs10 <- c(coeffs, 6) # 10次の項(係数6)を追加した、10次多項式の係数

f10 <- function(xi) {
  total <- 0
  for (k in 0:10) {
    total <- total + coeffs10[k + 1] * xi^k
  }
  total
}

analytic_integral10 <- 0
for (k in 0:10) {
  if (k %% 2 == 0) {
    analytic_integral10 <- analytic_integral10 + coeffs10[k + 1] * (2 / (k + 1))
  }
}

gauss_integral10 <- sum(gauss_wts * f10(gauss_pts))

cat("--- 10次多項式(5点ガウス求積の適用範囲を超える場合) ---\n")
cat("解析的な積分値:", analytic_integral10, "\n")
cat("5点ガウス求積の近似値:", gauss_integral10, "\n")
cat("差:", gauss_integral10 - analytic_integral10, "\n")