1. 問題設定（どんな方程式を解くのか）

ここでは，もっとも基本的な一次常微分方程式を扱います。

対象とする方程式は

\[ \frac{dx}{dt} = -kx \]

初期条件は

\[ x(0)=x_0 \]

です。

この式は，「現在の値 \(x\) に比例して減少する」ことを意味しており，
物理・化学・地球科学などで多く登場する 指数関数的減衰 を表します
（例：放射性崩壊，電荷の放電，冷却，減衰振動の外側など）。

今回は次の値を用いて比較します：

library(deSolve)
library(ggplot2)

x0    <- 5
k     <- 0.5
t_end <- 10
dt    <- 0.1
times <- seq(0, t_end, by = dt)

2. 解析解（手計算で求められる正確な解）

この方程式は変数分離によって解析的に解くことができます。

\[ \frac{dx}{dt}=-kx \quad \Rightarrow \quad x(t)=x_0 e^{-kt} \]

まず，この解析解を基準として計算します。

sol_analytic <- function(t) x0 * exp(-k * t)

exact <- data.frame(
  time = times,
  x    = sol_analytic(times),
  type = "解析解"
)

3. オイラー法（基本的な数値解法）

解析的に解けない微分方程式も多く存在するため，
数値的に解く方法が非常に重要になります。

最も基本的な方法が オイラー法 です。

微分方程式

\[ \frac{dx}{dt}=f(t,x) \]

に対して，離散化すると

\[ x_{n+1} = x_n + \Delta t\, f(t_n, x_n) \]

となります。

今回の場合

\[ f(t,x) = -kx \]

なので，
\[ x_{n+1}=x_n - k\Delta t x_n \]

です。

euler_decay <- function(x0, k, dt, t_end) {
  n <- floor(t_end / dt)
  t <- seq(0, n * dt, by = dt)
  x <- numeric(length(t))
  x[1] <- x0
  
  for (i in 1:n) {
    dx_dt  <- -k * x[i]
    x[i+1] <- x[i] + dt * dx_dt
  }
  
  data.frame(time = t, x = x, type = "オイラー法")
}

res_euler <- euler_decay(x0, k, dt, t_end)

4. deSolve による数値解（実用的な解法）

deSolve パッケージは，
微分方程式をより高精度で効率的に解くための R の代表的なツールです。

内部では Runge–Kutta 法 や 線形多段法 などが柔軟に選択されるため，
解析解が存在しない複雑な微分方程式でも非常に強力に機能します。

ここでは，オイラー法と比較する目的で ode() を使います。

test <- function(t, state, parms) {
  x <- state[1]
  k <- parms[1]
  dx <- -k * x
  list(c(dx))
}

initial <- c(x = x0)
parms   <- c(k)

res_ode <- ode(y = initial, times = times, func = test, parms = parms)
res_ode <- as.data.frame(res_ode)
res_ode$type <- "deSolve"

5. 3つの解を比較（解析解・オイラー法・deSolve）

3つの解を同じ図にまとめて比較します。

plot_all <- rbind(
  exact,
  res_euler,
  res_ode
)

ggplot(plot_all, aes(x = time, y = x, linetype = type, shape = type, col=type)) +
  geom_line(data = subset(plot_all, type == "解析解"), linewidth = 0.5) +
  geom_point(data = subset(plot_all, type != "解析解"), size = 1.5) +
  scale_colour_manual(values = c( "解析解"   = "black",  "オイラー法" = "red", "deSolve" = "blue")) +
  labs(
    title = "解析解 vs オイラー法 vs deSolve",
    x = "Time",
    y = "x(t)"
  ) +
  theme(
    panel.background = element_rect(fill="white", colour="black", linetype="solid"),
    plot.title = element_text(size=18, face="bold", hjust=0.5),
    axis.title = element_text(size=16),
    axis.text  = element_text(size=13),
    legend.text = element_text(size=13)
  )

6. 誤差評価（オイラー法と deSolve を同時比較）

解析解との差の絶対誤差

\[ | x_{\text{num}} - x_{\text{exact}} | \]

を計算し，オイラー法と deSolve を同じ図に重ねて比較します。

df_exact <- data.frame(time = exact$time, x_exact = exact$x)
df_euler <- data.frame(time = res_euler$time, x_euler = res_euler$x)
df_ode   <- data.frame(time = res_ode$time,   x_deSolve = res_ode$x)

merged_all <- merge(df_exact, df_euler, by="time")
merged_all <- merge(merged_all, df_ode,   by="time")

merged_all$abs_error_Euler   <- abs(merged_all$x_euler   - merged_all$x_exact)
merged_all$abs_error_deSolve <- abs(merged_all$x_deSolve - merged_all$x_exact)

err_long <- rbind(
  data.frame(time=merged_all$time, abs_error=merged_all$abs_error_Euler, method="オイラー法"),
  data.frame(time=merged_all$time, abs_error=merged_all$abs_error_deSolve, method="deSolve")
)

ggplot(err_long, aes(x=time, y=abs_error, col=method)) +
  geom_line(linewidth=1.2) +
  geom_point(size=3) +
  scale_colour_manual(values = c( "解析解"   = "black",  "オイラー法" = "red", "deSolve" = "blue")) +
  labs(
    title="オイラー法と deSolve の誤差",
    x="Time",
    y="絶対誤差 |numerical - exact|"
  ) +
  theme(
    panel.background = element_rect(fill="white", colour="black", linetype="solid"),
    plot.title = element_text(size=18, face="bold", hjust=0.5),
    axis.title = element_text(size=16),
    axis.text  = element_text(size=13),
    legend.text = element_text(size=13)
  )

7. まとめ