はじめに

　R言語を使って、円関数の定義や公式を可視化して理解しようシリーズです。

　この記事では、対数渦巻き(ベルヌーイの渦巻き)のグラフを作成します。

【前の内容】

www.anarchive-beta.com

【他の記事一覧】

www.anarchive-beta.com

【この記事の内容】

はじめに
対数螺旋の可視化
参考書籍
おわりに

対数螺旋の可視化

　sin関数(サイン関数・sine function)とcos関数(コサイン関数・cosine function)を用いて定義される対数螺旋(対数渦巻・logarithmic spiral・ベルヌーイの螺旋・Bernoulli's spiral)をグラフで確認します。

　利用するパッケージを読み込みます。

# 利用パッケージ
library(tidyverse)
library(gganimate)

　この記事では、基本的に パッケージ名::関数名() の記法を使うので、パッケージの読み込みは不要です。ただし、作図コードについてはパッケージ名を省略するので、 ggplot2 を読み込む必要があります。
　また、ネイティブパイプ演算子 |> を使います。magrittr パッケージのパイプ演算子 %>% に置き換えられますが、その場合は magrittr を読み込む必要があります。

定義式の確認

　まずは、対数螺旋の定義を数式で確認します。
　sin関数については「【R】sin関数の可視化 - からっぽのしょこ」、cos関数については「【R】cos関数の可視化 - からっぽのしょこ」、度数法と弧度法の角度については「【R】円周の作図 - からっぽのしょこ」を参照してください。

　対数螺旋は、sin関数とcos関数を用いて、次の式で定義されます。

$\displaystyle \begin{align} &\quad r = a e^{b \theta} \\ &\left\{ \begin{aligned} x &= r \cos \theta \\ y &= r \sin \theta \end{aligned} \right. \end{align}$

　変数 $\theta$ は弧度法における角度(ラジアン)です。度数法における角度を $\theta^{\circ}$ とすると、 $\theta = \frac{2 \pi}{360^{\circ}} \theta^{\circ}$ の関係です。 $e$ はネイピア数(自然対数の底)、 $\pi$ は円周率です。
　螺旋の形状は定数(パラメータ) $a, b$ によって決まります。正の実数 $a \gt 0$ と実数 $b$ は螺旋の間隔や全体のサイズに影響します。
　 $r$ は中心からの距離で、 $\theta$ に応じて変化します。

　 $r$ の式を $\theta$ について整理します。

$\displaystyle \begin{align} r &= a e^{b \theta} \\ \Rightarrow \log r &= \log a + b \theta \\ \Rightarrow b \theta &= \log r - \log a \\ &= \log \frac{r}{a} \\ \Rightarrow \theta &= \frac{1}{b} \log \frac{r}{a} \end{align}$

　両辺の対数をとり、対数の性質 $\log e^x = x$ 、 $\log (x y) = \log x + \log y$ 、 $\log \frac{x}{y} = \log x - \log y$ を用いて式変形しています。
　 $r$ の式は指数ですが、 $\theta$ の式は対数で表されます。

対数螺旋の作図

　次は、対数螺旋のグラフを作成して、変数(ラジアン)と座標(曲線上の点)の関係を確認します。

螺旋の作図

　対数螺旋のグラフを作成します。

　パラメータを指定して、螺旋の描画用のデータフレームを作成します。

# 定数を指定
a <- 1
b <- 0.1

# 周回数を指定
lap_num <- 4

# 螺旋用のラジアン(弧度法の角度)を作成
t_vec <- seq(from = 0, to = lap_num*2*pi, length.out = 1000) # 正の範囲
#t_vec <- seq(from = -lap_num*2*pi, to = 0, length.out = 1000) # 負の範囲
#t_vec <- seq(from = -lap_num*pi, to = lap_num*pi, length.out = 1000) # 正と負の範囲

# 螺旋の座標を作成
spiral_df <- tibble::tibble(
  t = t_vec, # ラジアン
  r = a * exp(b * t), # ノルム
  x = r * cos(t), # x座標
  y = r * sin(t), # y座標
  sgn_t_flag = t >= 0 # 角度の正負
)
spiral_df

## # A tibble: 1,000 × 5
##         t     r     x      y sgn_t_flag
##     <dbl> <dbl> <dbl>  <dbl> <lgl>     
##  1 0       1    1     0      TRUE      
##  2 0.0252  1.00 1.00  0.0252 TRUE      
##  3 0.0503  1.01 1.00  0.0505 TRUE      
##  4 0.0755  1.01 1.00  0.0760 TRUE      
##  5 0.101   1.01 1.01  0.101  TRUE      
##  6 0.126   1.01 1.00  0.127  TRUE      
##  7 0.151   1.02 1.00  0.153  TRUE      
##  8 0.176   1.02 1.00  0.178  TRUE      
##  9 0.201   1.02 1.00  0.204  TRUE      
## 10 0.226   1.02 0.997 0.230  TRUE      
## # ℹ 990 more rows

　周回数 lap_num を指定して、螺旋の座標計算用のラジアン $\theta$ を作成します。 $2 \pi$ の範囲で1周します。
　実数 $a$ を指定して、原点からの距離(変数ごとの係数) $r = a e^{b \theta}$ を計算し、x軸の値 $x = r \cos \theta$ とy軸の値 $y = r \sin \theta$ を計算します。

　グラフサイズや軸線用の値を設定します。

# グラフサイズを設定
axis_size <- c(spiral_df[["x"]], spiral_df[["y"]]) |> 
  abs() |> 
  max() |> 
  ceiling()

# 目盛間隔を設定
step_val <- 5

# 軸線数を設定
if(axis_size%%step_val == 0) {
  
  # グラフサイズと最大目盛が一致する場合
  circle_num <- axis_size %/% step_val
} else {
  
  # グラフサイズと最大目盛が一致しない場合
  circle_num <- axis_size %/% step_val + 1
  
  # グラフサイズを拡大
  axis_size <- max(axis_size, circle_num*step_val)
}
axis_size; circle_num

## [1] 15
## [1] 3

　ノルム軸線の目盛間隔(円ごとの半径の間隔)を step_val、軸線数(円の数)を circle_num とします。作図時に自動で設定されるx軸・y軸の目盛間隔と同じ値にすると見栄えがよくなります。
　螺旋のx軸・y軸の最大値とノルム軸の最大値からグラフサイズの半分を axis_size とします。

　ノルム軸線(円)の描画用のデータフレームを作成します。

# ノルム軸線の座標を作成
coord_circle_df <- tidyr::expand_grid(
  r = 1:circle_num * step_val, 
  t = seq(from = 0, to = 2*pi, length.out = 361), 
) |> # 半径ごとにラジアンを複製
  dplyr::mutate(
    x = r * cos(t), 
    y = r * sin(t)
  )
coord_circle_df

## # A tibble: 1,083 × 4
##        r      t     x      y
##    <dbl>  <dbl> <dbl>  <dbl>
##  1     5 0       5    0     
##  2     5 0.0175  5.00 0.0873
##  3     5 0.0349  5.00 0.174 
##  4     5 0.0524  4.99 0.262 
##  5     5 0.0698  4.99 0.349 
##  6     5 0.0873  4.98 0.436 
##  7     5 0.105   4.97 0.523 
##  8     5 0.122   4.96 0.609 
##  9     5 0.140   4.95 0.696 
## 10     5 0.157   4.94 0.782 
## # ℹ 1,073 more rows

　半径が step_val 間隔の circle_num 個の円周の座標を計算します。外側の円周の半径は axis_size になります。

　角度軸線(斜線)の描画用のデータフレームを作成します。

# 半円における目盛数(分母の値)を指定
denom <- 6

# 角度軸線の座標を作成
coord_oblique_df <- tibble::tibble(
  i = seq(from = 0, to = 2*denom-1, by = 1), # 目盛位置番号(分子の値)
  t = i / denom * pi, 
  r = axis_size, # ノルム軸線の最大値
  x = r * cos(t), 
  y = r * sin(t), 
  t_label = paste0("frac(", i, ", ", denom, ")~pi"), # 角度ラベル
  h = 1 - (x/r * 0.5 + 0.5), 
  v = 1 - (y/r * 0.5 + 0.5)
)
coord_oblique_df

## # A tibble: 12 × 8
##        i     t     r         x         y t_label             h      v
##    <dbl> <dbl> <dbl>     <dbl>     <dbl> <chr>           <dbl>  <dbl>
##  1     0 0        15  1.5 e+ 1  0        frac(0, 6)~pi  0      0.5   
##  2     1 0.524    15  1.30e+ 1  7.5 e+ 0 frac(1, 6)~pi  0.0670 0.25  
##  3     2 1.05     15  7.5 e+ 0  1.30e+ 1 frac(2, 6)~pi  0.25   0.0670
##  4     3 1.57     15  9.18e-16  1.5 e+ 1 frac(3, 6)~pi  0.5    0     
##  5     4 2.09     15 -7.5 e+ 0  1.30e+ 1 frac(4, 6)~pi  0.75   0.0670
##  6     5 2.62     15 -1.30e+ 1  7.5 e+ 0 frac(5, 6)~pi  0.933  0.25  
##  7     6 3.14     15 -1.5 e+ 1  1.84e-15 frac(6, 6)~pi  1      0.5   
##  8     7 3.67     15 -1.30e+ 1 -7.5 e+ 0 frac(7, 6)~pi  0.933  0.75  
##  9     8 4.19     15 -7.50e+ 0 -1.30e+ 1 frac(8, 6)~pi  0.75   0.933 
## 10     9 4.71     15 -2.76e-15 -1.5 e+ 1 frac(9, 6)~pi  0.5    1     
## 11    10 5.24     15  7.5 e+ 0 -1.30e+ 1 frac(10, 6)~pi 0.25   0.933 
## 12    11 5.76     15  1.30e+ 1 -7.50e+ 0 frac(11, 6)~pi 0.0670 0.75

　ノルム軸線の外側の円周上を等間隔の 2 * denom 個に分割した点の座標を計算します。denom を $n$ とすると、目盛間隔は $\frac{1}{n} \pi$ になります。また $i, n$ を整数として、 $\frac{i}{n} \pi$ と $\frac{i + 2 n}{n} \pi = \frac{i}{n} \pi + 2 \pi$ の目盛位置が一致します。

　螺旋のグラフを作成します。

# ラベル用の文字列を作成
fnc_label <- paste0(
  "list(", 
  "r == a * e^{b * theta}", ", ", 
  "a == ", a, ", ", 
  "b == ", b, 
  ")"
)

# 螺旋を作図
add_size <- 0.5
ggplot() + 
  geom_path(data = coord_circle_df, 
            mapping = aes(x = x, y = y, group = r), 
            color = "white") + # ノルム軸線
  geom_segment(data = coord_oblique_df, 
               mapping = aes(x = 0, y = 0, xend = x, yend = y, group = i), 
               color = "white") + # 角度軸線
  geom_text(data = coord_oblique_df, 
            mapping = aes(x = x, y = y, label = t_label, hjust = h, vjust = v), 
            parse = TRUE) + # 角度目盛ラベル
  geom_path(data = spiral_df, 
            mapping = aes(x = x, y = y), 
            linewidth = 1) + # 螺旋
  coord_fixed(ratio = 1, 
              xlim = c(-axis_size-add_size, axis_size+add_size), 
              ylim = c(-axis_size-add_size, axis_size+add_size)) + 
  labs(title = "logarithmic spiral", 
       subtitle = parse(text = fnc_label), 
       x = expression(x == r ~ cos~theta), 
       y = expression(y == r ~ sin~theta))

　(x軸方向に値が増減するため geom_line() ではなく、) geom_path() で螺旋を描画します。

　以上で、対数螺旋を描画できました。次からは、変数やパラメータによる螺旋への影響を確認していきます。

角度と座標の関係

　対数螺旋上を移動する点のアニメーションを作成します。

・作図コード(クリックで展開)

　フレーム数を指定して、螺旋上の点の描画用のデータフレームを作成します。

# フレーム数を指定
frame_num <- 300

# 螺旋上の点用のラジアンを作成
t_vals <- seq(from = min(t_vec), to = max(t_vec), length.out = frame_num+1)[1:frame_num]

# 螺旋上の点の座標を作成
anim_point_df <- tibble::tibble(
  frame_i = 1:frame_num, # フレーム番号
  t = t_vals, 
  r = a * exp(b * t), 
  x = r * cos(t), 
  y = r * sin(t), 
  var_label = paste0(
    "list(", 
    "a == ", a, ", ", 
    "b == ", b, ", ", 
    "theta == ", round(t/pi, digits = 2), " * pi, ", 
    "r == ", round(r, digits = 2), 
    ")"
  ) # 変数ラベル
)
anim_point_df

## # A tibble: 300 × 6
##    frame_i      t     r     x      y var_label                                  
##      <int>  <dbl> <dbl> <dbl>  <dbl> <chr>                                      
##  1       1 0       1    1     0      list(a == 1, b == 0.1, theta == 0 * pi, r …
##  2       2 0.0838  1.01 1.00  0.0844 list(a == 1, b == 0.1, theta == 0.03 * pi,…
##  3       3 0.168   1.02 1.00  0.170  list(a == 1, b == 0.1, theta == 0.05 * pi,…
##  4       4 0.251   1.03 0.993 0.255  list(a == 1, b == 0.1, theta == 0.08 * pi,…
##  5       5 0.335   1.03 0.977 0.340  list(a == 1, b == 0.1, theta == 0.11 * pi,…
##  6       6 0.419   1.04 0.953 0.424  list(a == 1, b == 0.1, theta == 0.13 * pi,…
##  7       7 0.503   1.05 0.921 0.507  list(a == 1, b == 0.1, theta == 0.16 * pi,…
##  8       8 0.586   1.06 0.883 0.587  list(a == 1, b == 0.1, theta == 0.19 * pi,…
##  9       9 0.670   1.07 0.838 0.664  list(a == 1, b == 0.1, theta == 0.21 * pi,…
## 10      10 0.754   1.08 0.786 0.738  list(a == 1, b == 0.1, theta == 0.24 * pi,…
## # ℹ 290 more rows

　フレーム数 frame_num を指定して、螺旋の座標計算用のラジアン t_vec の範囲を等間隔に frame_num 個に分割して、螺旋上の点の座標計算用のラジアン t_vals を作成します。

　角度を示す線分の描画用のデータフレームを作成します。

# 角度線の座標を作成
anim_angle_oblique_df <- tibble::tibble(
  frame_i = 1:frame_num, 
  t = t_vals, 
  r = a * exp(b * t), 
  x = axis_size * cos(t), 
  y = axis_size * sin(t)
)
anim_angle_oblique_df

## # A tibble: 300 × 5
##    frame_i      t     r     x     y
##      <int>  <dbl> <dbl> <dbl> <dbl>
##  1       1 0       1     15    0   
##  2       2 0.0838  1.01  14.9  1.26
##  3       3 0.168   1.02  14.8  2.50
##  4       4 0.251   1.03  14.5  3.73
##  5       5 0.335   1.03  14.2  4.93
##  6       6 0.419   1.04  13.7  6.10
##  7       7 0.503   1.05  13.1  7.23
##  8       8 0.586   1.06  12.5  8.30
##  9       9 0.670   1.07  11.8  9.32
## 10      10 0.754   1.08  10.9 10.3 
## # ℹ 290 more rows

　 $\theta$ に応じて、半径が axis_size の円周上の点の座標を計算します。

　螺旋と角度線上の点の描画用のデータフレームを作成します。

# 角度線上の点の座標を作成
r_min <- a * exp(b * min(t_vec))
r_max <- a * exp(b * max(t_vec))
anim_angle_point_df <- tidyr::expand_grid(
  frame_i = 1:frame_num, 
  lap_i   = (-lap_num):(lap_num-1) # 点番号
) |> # フレームごとに点を複製
  dplyr::mutate(
    tmp_t = t_vals[frame_i] %% (2*pi), # 単位円相当のラジアン
    t = lap_i * 2*pi + tmp_t, 
    r = a * exp(b * t), 
    x = r * cos(t), 
    y = r * sin(t),   
    sgn_t_flag = t >= 0 # 角度の正負
  ) |> 
  dplyr::filter(r >= r_min, r <= r_max) # 螺旋上の点を抽出
anim_angle_point_df

## # A tibble: 1,200 × 8
##    frame_i lap_i  tmp_t       t     r     x         y sgn_t_flag
##      <int> <int>  <dbl>   <dbl> <dbl> <dbl>     <dbl> <lgl>     
##  1       1     0 0       0       1     1     0        TRUE      
##  2       1     1 0       6.28    1.87  1.87 -4.59e-16 TRUE      
##  3       1     2 0      12.6     3.51  3.51 -1.72e-15 TRUE      
##  4       1     3 0      18.8     6.59  6.59 -4.84e-15 TRUE      
##  5       2     0 0.0838  0.0838  1.01  1.00  8.44e- 2 TRUE      
##  6       2     1 0.0838  6.37    1.89  1.88  1.58e- 1 TRUE      
##  7       2     2 0.0838 12.7     3.54  3.53  2.96e- 1 TRUE      
##  8       2     3 0.0838 18.9     6.64  6.62  5.56e- 1 TRUE      
##  9       3     0 0.168   0.168   1.02  1.00  1.70e- 1 TRUE      
## 10       3     1 0.168   6.45    1.91  1.88  3.18e- 1 TRUE      
## # ℹ 1,190 more rows

　 $\theta$ の前後 $2 \pi$ 間隔 $\theta + 2 i \pi$ の点の座標を計算します。 $i$ は整数で、lap_i 列とします。
　lap_i 列は、 $r$ が正の値のみの場合は 0 から lap_num - 1、負の値のみの場合は -lap_num から -1、正負の値を含む場合は ceiling(-0.5 * lap_num) から floor(0.5 * lap_num - 1) の整数です。(全ての場合を共通のコードで処理するために、)-lap_num から lap_num - 1 の整数を用いて座標を計算し、r の最小値から最大値の点(行)を取り出します。

　螺旋上の点のアニメーションを作成します。

# 螺旋上の点のアニメーションを作図
add_size <- 0.5
anim <- ggplot() + 
  geom_path(data = coord_circle_df, 
            mapping = aes(x = x, y = y, group = r), 
            color = "white") + # ノルム軸線
  geom_segment(data = coord_oblique_df, 
               mapping = aes(x = 0, y = 0, xend = x, yend = y, group = i), 
               color = "white") + # 角度軸線
  geom_segment(mapping = aes(x = c(-Inf, 0), y = c(0, -Inf), 
                             xend = c(Inf, 0), yend = c(0, Inf)), 
               arrow = arrow(length = unit(10, units = "pt"), ends = "last")) + # x・y軸線
  geom_text(data = coord_oblique_df, 
            mapping = aes(x = x, y = y, label = t_label, hjust = h, vjust = v), 
            parse = TRUE) + # 角度目盛ラベル
  geom_path(data = spiral_df, 
            mapping = aes(x = x, y = y, color = sgn_t_flag), 
            linewidth = 0.5) + # 螺旋
  geom_segment(data = anim_angle_oblique_df,
               mapping = aes(x = 0, y = 0, xend = x, yend = y),
               linewidth = 1) + # 角度用の補助線
  geom_point(data = anim_point_df, 
             mapping = aes(x = x, y = y), 
             size = 5) + # 螺旋上の点
  geom_point(data = anim_angle_point_df, 
             mapping = aes(x = x, y = y, color = sgn_t_flag), 
             size = 4, shape = "circle open") + # 角度用の補助線上の点
  geom_text(data = anim_point_df, 
            mapping = aes(x = -Inf, y = Inf, label = var_label), 
            parse = TRUE, hjust = 0, vjust = -0.5) + # 変数ラベル
  gganimate::transition_manual(frames = frame_i) + # フレーム切替
  scale_color_manual(breaks = c("TRUE", "FALSE"), 
                     values = c("blue", "red"), 
                     labels = c(expression(theta >= 0), expression(theta < 0)), 
                     name = expression(sgn~theta)) + # 角度の正負
  coord_fixed(ratio = 1, clip = "off", 
              xlim = c(-axis_size-add_size, axis_size+add_size), 
              ylim = c(-axis_size-add_size, axis_size+add_size)) + 
  labs(title = "logarithmic spiral", 
       subtitle = "", # (変数ラベルの表示用)
       x = expression(x == r ~ cos~theta), 
       y = expression(y == r ~ sin~theta))

# gif画像を作成
gganimate::animate(
  plot = anim, 
  nframes = frame_num, fps = 30, 
  width = 600, height = 600, 
  renderer = gganimate::gifski_renderer()
)

　gganimate パッケージを利用してアニメーション(gif画像)を作成します。
　transition_manual() のフレーム制御の引数 frames にフレーム番号列 frame_i を指定します。
　animate() の plot 引数にグラフオブジェクト、nframes 引数にフレーム数 frame_num を指定して、gif画像を作成します。また、fps 引数に1秒当たりのフレーム数を指定できます。

　gganimate を利用したフレーム切り替えでは、フレームに応じてサブタイトルを変更できますが、その際に expression() を使えないので、geom_text() を使って疑似的にサブタイトルの位置にパラメータラベルを表示しています。
　coord_***() に clip = "off" を指定すると、描画領域外(余白領域)に描画できます。

　フレームごとの $\theta$ の点を黒色の点、原点と $\theta$ の点を通る半直線を黒色の線分で示します。
　x軸線の正の部分と半直線の偏角が変数 $\theta$ に対応します。
　 $\theta + 2 i \pi$ の点( $i$ を整数として $\theta$ の前後 $2 \pi$ 間隔の点)が線分上に並ぶのを確認できます。

パラメータと形状の関係

　続いて、定数(パラメータ)の値の変化に応じたグラフを作成して、パラメータと螺旋の形状の関係を確認します。

・作図コード(クリックで展開)

　フレーム数とパラメータを指定して、螺旋の描画用のデータフレームを作成します。

# フレーム数を指定
frame_num <- 101

# パラメータを指定
a_vals <- seq(from = 0, to = 5, length.out = frame_num)
b_vals <- seq(from = 0.1, to = 0.1, length.out = frame_num)

# 周回数を指定
lap_num <- 10

# 点数を指定
point_num <- lap_num * 12 + 1

# 螺旋の座標を作成
anim_spiral_df <- tidyr::expand_grid(
  frame_i = 1:frame_num, 
  t       = seq(from = 0, to = lap_num*2*pi, length.out = 1000)
) |> # フレームごとにラジアンを複製
  dplyr::mutate(
    a = a_vals[frame_i], 
    b = b_vals[frame_i], 
    r = a * exp(b * t), 
    x = r * cos(t), 
    y = r * sin(t)
  )
anim_spiral_df

## # A tibble: 101,000 × 7
##    frame_i      t     a     b     r     x     y
##      <int>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1       1 0          0   0.1     0     0     0
##  2       1 0.0629     0   0.1     0     0     0
##  3       1 0.126      0   0.1     0     0     0
##  4       1 0.189      0   0.1     0     0     0
##  5       1 0.252      0   0.1     0     0     0
##  6       1 0.314      0   0.1     0     0     0
##  7       1 0.377      0   0.1     0     0     0
##  8       1 0.440      0   0.1     0     0     0
##  9       1 0.503      0   0.1     0     0     0
## 10       1 0.566      0   0.1     0     0     0
## # ℹ 100,990 more rows

　フレーム数 frame_num を指定して、frame_num 個の値を指定します。
　フレーム番号( 1 から frame_num までの整数)とラジアンの値の組み合わせを expand_grid() で作成します。パラメータの値ごとに螺旋用のラジアンを複製できます。

　目安として用いる螺旋上の点の描画用のデータフレームを作成します。

# 螺旋上の点の座標を作成
anim_point_df <- tidyr::expand_grid(
  frame_i = 1:frame_num, 
  t       = seq(
    from = min(anim_spiral_df[["t"]]), 
    to   = max(anim_spiral_df[["t"]]), 
    length.out = point_num
  )
) |> # フレームごとにラジアンを複製
  dplyr::mutate(
    a = a_vals[frame_i], 
    b = b_vals[frame_i], 
    r = a * exp(b * t), 
    x = r * cos(t), 
    y = r * sin(t)
  )
anim_point_df

## # A tibble: 12,221 × 7
##    frame_i     t     a     b     r     x     y
##      <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1       1 0         0   0.1     0     0     0
##  2       1 0.524     0   0.1     0     0     0
##  3       1 1.05      0   0.1     0     0     0
##  4       1 1.57      0   0.1     0     0     0
##  5       1 2.09      0   0.1     0     0     0
##  6       1 2.62      0   0.1     0     0     0
##  7       1 3.14      0   0.1     0     0     0
##  8       1 3.67      0   0.1     0     0     0
##  9       1 4.19      0   0.1     0     0     0
## 10       1 4.71      0   0.1     0     0     0
## # ℹ 12,211 more rows

　螺旋の対応関係の確認用に、変数(ラジアン)に関して一定間隔に点を描画することにします。
　点の数 point_num を指定し point_num 個のラジアンの値を作成して、螺旋の座標計算と同様に処理します。

　パラメータラベルの描画用のデータフレームを作成します。

# パラメータラベル用の文字列を作成
anim_label_df <- tibble::tibble(
  frame_i = 1:frame_num, # フレーム番号
  a       = a_vals, 
  b       = b_vals, 
  param_label = paste0(
    "list(", 
    "a == ", round(a, digits = 2), ", ", 
    "b == ", round(b, digits = 2), 
    ")"
  ) # パラメータラベル
)
anim_label_df

## # A tibble: 101 × 4
##    frame_i     a     b param_label              
##      <int> <dbl> <dbl> <chr>                    
##  1       1  0      0.1 list(a == 0, b == 0.1)   
##  2       2  0.05   0.1 list(a == 0.05, b == 0.1)
##  3       3  0.1    0.1 list(a == 0.1, b == 0.1) 
##  4       4  0.15   0.1 list(a == 0.15, b == 0.1)
##  5       5  0.2    0.1 list(a == 0.2, b == 0.1) 
##  6       6  0.25   0.1 list(a == 0.25, b == 0.1)
##  7       7  0.3    0.1 list(a == 0.3, b == 0.1) 
##  8       8  0.35   0.1 list(a == 0.35, b == 0.1)
##  9       9  0.4    0.1 list(a == 0.4, b == 0.1) 
## 10      10  0.45   0.1 list(a == 0.45, b == 0.1)
## # ℹ 91 more rows

　フレームごとにパラメータラベル用の文字列を作成します。

　グラフサイズや軸線用の値を設定します。

# グラフサイズを指定
axis_size <- 100

# 目盛間隔を指定
step_val <- 25

# 軸線数を設定
circle_num <- (axis_size * sqrt(2)) %/% step_val
circle_num

## [1] 5

　グラフサイズ(の半分の値) axis_size、ノルム軸線の目盛間隔 step_val を指定して、軸線数 circle_num を計算します。軸線数は、対角方向(縦横ではなく斜め)のグラフサイズから求めます。
　この例では、x軸・y軸とノルム軸の最大値を指定して、螺旋などが超える場合は表示しない(枠外になる)ことにします。

　ノルム軸線と角度軸線の描画用のデータフレームを作成します。

# ノルム軸線の座標を作成
coord_circle_df <- tidyr::expand_grid(
  r = 1:circle_num * step_val, 
  t = seq(from = 0, to = 2*pi, length.out = 361), 
) |> # 半径ごとにラジアンを複製
  dplyr::mutate(
    x = r * cos(t), 
    y = r * sin(t)
  )

# 半円における目盛数を指定
oblique_num <- 6

# 角度軸線の座標を作成
coord_oblique_df <- tibble::tibble(
  i = 0:(oblique_num-1), # 目盛位置番号
  t = i / oblique_num * pi, 
  a = tan(t) # ノルム軸線の傾き
)
coord_oblique_df

## # A tibble: 6 × 3
##       i     t         a
##   <int> <dbl>     <dbl>
## 1     0 0      0       
## 2     1 0.524  5.77e- 1
## 3     2 1.05   1.73e+ 0
## 4     3 1.57   1.63e+16
## 5     4 2.09  -1.73e+ 0
## 6     5 2.62  -5.77e- 1

　ノルム軸線用の値は、「螺旋の作図」のときと同様にして作成します。ただし、ノルム軸の最大値の設定の仕方が変わっています。
　角度軸線用の値は、 $0 \leq \theta \lt \pi$ の範囲における軸線数(1周における軸線数の半数)を $n$ 、軸線番号を $i$ として、軸線の傾き $a = \tan (\frac{i}{n} \pi)$ を用います。

　螺旋のアニメーションを作成します。

# ラベル用の文字列を作成
fnc_label <- "r == a * e^{b * theta}"

# 螺旋のアニメーションを作図
anim <- ggplot() + 
  geom_path(data = coord_circle_df, 
            mapping = aes(x = x, y = y, group = r), 
            color = "white") + # ノルム軸線
  geom_abline(data = coord_oblique_df, 
              mapping = aes(slope = a, intercept = 0), 
              color = "white") + # 角度軸線
  geom_path(data = anim_spiral_df, 
            mapping = aes(x = x, y = y, color = t/pi), 
            linewidth = 1) + # 螺旋
  geom_point(data = anim_point_df, 
             mapping = aes(x = x, y = y, color = t/pi), 
             size = 2) + # 螺旋上の点
  geom_label(data = anim_label_df, 
             mapping = aes(x = -Inf, y = Inf, label = param_label), 
             parse = TRUE, hjust = 0, vjust = 1, alpha = 0.5, label.r = unit(0, units = "pt")) + # パラメータラベル
  gganimate::transition_manual(frames = frame_i) + # フレーム切替
  coord_fixed(ratio = 1, 
              xlim = c(-axis_size, axis_size), 
              ylim = c(-axis_size, axis_size)) + 
  labs(title = "logarithmic spiral", 
       subtitle = parse(text = fnc_label), 
       color = expression(frac(theta, pi)), 
       x = expression(x == r ~ cos~theta), 
       y = expression(y == r ~ sin~theta))

# gif画像を作成
gganimate::animate(
  plot = anim, 
  nframes = frame_num, fps = 10, 
  width = 600, height = 600, 
  renderer = gganimate::gifski_renderer()
)

　 $a$ の絶対値が大きいほど、 $r$ の絶対値が大きくなるため、螺旋のサイズが大きくなります。 $a$ の符号によってx軸・y軸に対して反転します。

対数螺旋とsin曲線・cos曲線の関係

　最後は、対数螺旋とsin関数・cos関数の曲線のグラフを作成して、変数(ラジアン)と座標や定数(パラメータ)と形状の関係を確認します。
　作図コードについては「GitHub - anemptyarchive/Mathematics: 数学ノート」、作図の解説については「【R】cos関数の可視化 - からっぽのしょこ」などを参照してください。

変数と座標の関係

　変数に応じて移動する対数螺旋・sin曲線・cos曲線上の点のアニメーションを作成します。

　螺旋上の点(黒色の点)のx軸方向の変化(青色の点)とy軸方向の変化(赤色の点)が、それぞれcos関数曲線上の点とsin関数曲線上の点に対応しているのを確認できます。
　 $\theta$ の絶対値が大きいほど、 $r$ の絶対値が大きくなるため、曲線の振幅が大きくなります。

パラメータと形状の関係

　パラメータに応じて変化する対数螺旋・sin曲線・cos曲線のアニメーションを作成します。

　対応関係の目安として、 $\theta$ に応じて曲線を色付けし、また一定間隔で点を打っています。螺旋のx軸方向の形状とy軸方向の形状が、それぞれcos関数曲線の形状とsin関数曲線の形状に対応しているのを確認できます。
　 $a$ の絶対値が大きいほど、 $r$ の絶対値が大きくなるため、曲線の振幅が大きくなります。 $a$ の符号によってそれぞれの軸に対して反転します。

　この記事では、対数螺旋のグラフを作成しました。次の記事では、フェルマーの螺旋のグラフを作成します。

参考書籍

『曲線と曲面（改訂版）－微分幾何的アプローチ－』梅原雅顕・山田光太郎,裳華房,2015年.

おわりに

　個人的に、この曲線が一番渦巻きって感じがします。対数じゃなくて指数螺旋ではと言いたくなるのは私だけではないようです。一応、式変形して対数に登場してもらいました。最近は数式弄りを避けていたので、久しぶりに式変形した気がします。
　円関数の文脈で円周率と共にiが出てくると虚数単位と思ってしまいますが、インデックスです。kは他の記事で係数として使ったので、jにした方がよかったかもしれません。円関数シリーズを書き始めた段階では、iを見て虚数単位と思うことはなかったので、これも成長なんでしょう。

　シリーズで記事を書く時は、章単位でまとめて書き終わってから1つずつ投稿していく方針にしたのですが、諸事情で1つ書いてはアップしています。この方法だと、書き進めるごとに書き換えたり書き足したりで、生産性の少ない修正が増えるのでもう止めたいのですが。今回もご多分に漏れず、もう1つ目の記事に修正したいことがあります、、例えば、円と斜線による追加の軸線の描き方とかです。
　そんな面倒事よりも優先されるって一体どんな事情なんだ。

　2023年10月20日は、BEYOOOOONDSの岡村美波さんの19歳のお誕生日です。

　みいみはただでさえ可愛いというのに、おそらく意識的にアイドルをされていてとっても可愛い方です。

【次の内容】

つづく

からっぽのしょこ

読んだら書く！書いたら読む！同じ事は二度調べ(たく)ない

【R】対数螺旋の可視化

はじめに

対数螺旋の可視化

定義式の確認

対数螺旋の作図

螺旋の作図

角度と座標の関係

パラメータと形状の関係

対数螺旋とsin曲線・cos曲線の関係

変数と座標の関係

パラメータと形状の関係

参考書籍

おわりに