LoginSignup
1
1

More than 1 year has passed since last update.

@zhidao(Zhidao)

一般化Fresnel積分の数値解法

Last updated at Posted at 2021-09-28

Fresnel積分の基本形

Fresnel（フレネル）積分とは
$$
\int_0^{x} e^{i\frac{\pi}{2}\sigma^{2}}\mathrm{d}\sigma=c(x)+is(x)
$$
または
$$
\begin{align}
\int_0^{x}\cos\frac{\pi}{2}\sigma^{2}\mathrm{d}\sigma&=c(x) \\
\int_0^{x}\sin\frac{\pi}{2}\sigma^{2}\mathrm{d}\sigma&=s(x)
\end{align}
$$
と書いた時の$c(x)$と$s(x)$の組のことです。これが何を意味するかについてはWikipediaに委ねます。いわゆる特殊関数に含まれているもので、解析的に計算することができません。ただし数値計算法はよく知られており、Numerical Recipesにも書かれています。

参考までに、ZMの中での実装コードを示します。

#define Z_FRESNEL_FP_MIN ( 1.0e-30 )

/* Fresnel integral by Gilbert formula */
static bool _zFresnelIntgGilbert(double x, double *s, double *c)
{
  int i, n;
  double sum, t, t0;

  *s = sum = 0;
  *c = t = x;
  t0 = zPI_2 * x * x;
  for( n=3, i=1; i<=Z_MAX_ITER_NUM; i++, n+=2 ){
    t *= t0/i;
    sum += ( i % 4 <= 1 ? t : -t ) / n;
    if( t < fabs(sum) * zTOL ) break;
    if( i % 2 ){
      *s = sum;
      sum = *c;
    } else{
      *c = sum;
      sum = *s;
    }
  }
  if( i > Z_MAX_ITER_NUM ){
    ZITERWARN( Z_MAX_ITER_NUM );
    return false;
  }
  return true;
}

/* Fresnel integral by modified Lentz method */
static bool _zFresnelIntgLentz(double x, double *s, double *c)
{
  int i, n;
  double a, pix2;
  zComplex z1, b, ca, cd, cc, d, h, del, cs;
  bool ret = true;

  pix2 = zPI_2 * x * x;
  _zComplexCreate( &z1, 1.0, 0.0 );
  _zComplexCreate( &b, 1.0, -pix2*2 );
  _zComplexCreate( &cc, 1.0/Z_FRESNEL_FP_MIN, 0.0 );
  zComplexCDiv( &z1, &b, &d );
  _zComplexCopy( &d, &h );
  for( n=1, i=2; i<=Z_MAX_ITER_NUM; i++, n+=2 ){
    a = -n * (n+1);
    b.re += 4.0;
    _zComplexMul( &d, a, &cd );
    _zComplexAddDRC( &cd, &b );
    zComplexCDiv( &z1, &cd, &d );
    zComplexInv( &cc, &cd );
    _zComplexMulDRC( &cd, a );
    _zComplexAdd( &cd, &b, &cc );
    _zComplexCMul( &cc, &d, &del );
    zComplexCMulDRC( &h, &del );
    del.re -= 1.0;
    if( zComplexIsTiny( &del ) ) break;
  }
  if( i > Z_MAX_ITER_NUM ){
    ZITERWARN( Z_MAX_ITER_NUM );
    ret = false;
  }
  _zComplexCreateCMul( &ca, h.re, h.im, x, -x );
  _zComplexCreateCMul( &cd, cos(pix2), sin(pix2), ca.re, ca.im );
  _zComplexCreateCMul( &cs, 0.5, 0.5, 1-cd.re, -cd.im );
  *c = cs.re;
  *s = cs.im;
  return ret;
}

/* Fresnel integral */
bool zFresnelIntgPI_2(double x, double *s, double *c)
{
  double ax;
  const double xth = 1.5; /* a heuristic threshould */
  bool ret = true;

  if( ( ax = fabs( x ) ) < sqrt( Z_FRESNEL_FP_MIN ) ){ /* avoid failure of convergence */
    *s = 0;
    *c = ax;
  } else if( ax <= xth ){ /* Gilbert formula for a small argument */
    ret = _zFresnelIntgGilbert( ax, s, c );
  } else{ /* modified Lentz method for a large argument */
    ret = _zFresnelIntgLentz( ax, s, c );
  }
  if( x < 0 ){ /* antisymmetry */
    *c = -*c;
    *s = -*s;
  }
  return ret;
}

複素数使った修正Lentz法の所がちょっと見づらいですが、ご勘弁を。

Fresnel積分の一般化

さて、上記のように既に知られている数値計算法を使い、次の計算をしたいというのがこの記事の本題です。
$$
\int_0^{x} \exp {i(\phi_{0}+\phi_{1}\sigma+\phi_{2}\sigma^{2})}\mathrm{d}\sigma=C(x)+iS(x)
$$
こんなもの何に使うの？という疑問は尤もで、クロソイド補間に使います。例えば次の牧野先生の方法とか。

基本形のスケーリング

まずは次を計算する方法から考えましょう。
$$
\int_0^{x} e^{i\phi_{2}\sigma^{2}}\mathrm{d}\sigma=C_{2}(x;\phi_{2})+iS_{2}(x;\phi_{2})
$$
混乱を避けるために、いったん元の基本形の媒介変数を$\sigma^{\prime}$、変数を$x^{\prime}$とそれぞれ置き直します。
$$
\int_0^{x^{\prime}} e^{i\frac{\pi}{2}\sigma^{\prime 2}}\mathrm{d}\sigma^{\prime}=c(x^{\prime})+is(x^{\prime})
$$
意味は変わっていません。$\phi_{2}$は任意の実数を取り得ますが、簡単のためここでは正である$(\phi_{2}>0)$と仮定しましょう。基本形を上記に揃えるため
$$
\phi_{2}\sigma^{2}=\frac{\pi}{2}\sigma^{\prime 2}
$$
なる変数変換を行えば、
$$
\sigma^{\prime}=\sqrt{\frac{2\phi_{2}}{\pi}}\sigma,\qquad
\mathrm{d}\sigma^{\prime}=\sqrt{\frac{2\phi_{2}}{\pi}}\mathrm{d}\sigma
$$
であり、また積分範囲は
$$
\begin{align}
\sigma^{\prime}&:0\rightarrow x^{\prime}
\\
\sigma&:0\rightarrow \sqrt{\frac{2\phi_{2}}{\pi}}x^{\prime}
\end{align}
$$
と対応しますので、
$$
\sqrt{\frac{2\phi_{2}}{\pi}}\int_0^{\sqrt{\frac{2\phi_{2}}{\pi}}x^{\prime}} e^{i\phi_{2}\sigma^{2}}\mathrm{d}\sigma=c(x^{\prime})+is(x^{\prime})
$$
が言えます。よって、
$$
\begin{align}
C_{2}(x;\phi_{2})&=\sqrt{\frac{\pi}{2\phi_{2}}}c\left(\sqrt{\frac{2\phi_{2}}{\pi}}x\right)
\\
S_{2}(x;\phi_{2})&=\sqrt{\frac{\pi}{2\phi_{2}}}s\left(\sqrt{\frac{2\phi_{2}}{\pi}}x\right)
\end{align}
$$
となります。

$\phi_{2}$が負の場合は、
$$
\begin{align}
C_{2}(x;\phi_{2})=\int_0^{x}\cos\phi_{2}\sigma^{2}\mathrm{d}\sigma&=\int_0^{x}\cos\left(-\phi_{2}\sigma^{2}\right)\mathrm{d}\sigma=C_{2}(x;-\phi_{2}) \\
S_{2}(x;\phi_{2})=\int_0^{x}\sin\phi_{2}\sigma^{2}\mathrm{d}\sigma&=-\int_0^{x}\sin\left(-\phi_{2}\sigma^{2}\right)\mathrm{d}\sigma=-S_{2}(x;-\phi_{2})
\end{align}
$$
とすれば良いです。

ZMでの実装は次の通りです。

/* a scaled version of Fresnel integral */
bool zFresnelIntgScale(double x, double f, double *s, double *c)
{
  bool ret;
  double r;

  if( zIsTiny( f ) ){ /* sin 0x = 0, cos 0x = 1 */
    *s = 0;
    *c = x;
    return true;
  }
  r = sqrt( zPI_2/fabs(f) );
  ret = zFresnelIntgPI_2( x/r, s, c );
  *s *= f > 0 ? r : -r;
  *c *= r;
  return ret;
}

一般化された形の平方完成

$\phi_{2}\neq 0$ならば、
$$
\begin{align}
\int_0^{x} \exp {i(\phi_{0}+\phi_{1}\sigma+\phi_{2}\sigma^{2})}\mathrm{d}\sigma
&=\int_0^{x} \exp i\left\{\phi_{2}\left(\sigma+\frac{\phi_{1}}{2\phi_{2}}\right)^{2}+\phi_{0}-\frac{\phi_{1}^{2}}{4\phi_{2}}\right\}\mathrm{d}\sigma
\\
&=\exp i\left(\phi_{0}-\frac{\phi_{1}^{2}}{4\phi_{2}}\right)
\int_0^{x} \exp i\phi_{2}\left(\sigma+\frac{\phi_{1}}{2\phi_{2}}\right)^{2}\mathrm{d}\sigma
\\
&=\exp i\left(\phi_{0}-\frac{\phi_{1}^{2}}{4\phi_{2}}\right)
\int_{\frac{\phi_{1}}{2\phi_{2}}}^{x+\frac{\phi_{1}}{2\phi_{2}}} e^{i\phi_{2}\sigma^{2}}\mathrm{d}\sigma
\\
&=\exp i\left(\phi_{0}-\frac{\phi_{1}^{2}}{4\phi_{2}}\right)\left(
\int_{0}^{x+\frac{\phi_{1}}{2\phi_{2}}} e^{i\phi_{2}\sigma^{2}}\mathrm{d}\sigma
-\int_{0}^{\frac{\phi_{1}}{2\phi_{2}}} e^{i\phi_{2}\sigma^{2}}\mathrm{d}\sigma
\right)
\end{align}
$$
なので、
$$
\begin{align}
\bar{c}&=\cos\left(\phi_{0}-\frac{\phi_{1}^{2}}{4\phi_{2}}\right) \\
\bar{s}&=\sin\left(\phi_{0}-\frac{\phi_{1}^{2}}{4\phi_{2}}\right) \\
\bar{C}(x)&=C_{2}\left(x+\frac{\phi_{1}}{2\phi_{2}};\phi_{2}\right)-C_{2}\left(\frac{\phi_{1}}{2\phi_{2}};\phi_{2}\right) \\
\bar{S}(x)&=S_{2}\left(x+\frac{\phi_{1}}{2\phi_{2}};\phi_{2}\right)-S_{2}\left(\frac{\phi_{1}}{2\phi_{2}};\phi_{2}\right)
\end{align}
$$
とおけば、
$$
C(x)+iS(x)=(\bar{c}+i\bar{s})(\bar{C}(x)+i\bar{S}(x))
$$
すなわち
$$
\begin{align}
C(x)&=\bar{c}\bar{C}(x)-\bar{s}\bar{S}(x) \\
S(x)&=\bar{c}\bar{S}(x)+\bar{s}\bar{C}(x)
\end{align}
$$
となります。

例外処理

$\phi_{2}=0$のときは、
$\phi_{1}\neq 0$なら
$$
\int_0^{x} \exp {i(\phi_{0}+\phi_{1}\sigma)}\mathrm{d}\sigma
=\frac{e^{i(\phi_{0}+\phi_{1}x)}-e^{i\phi_{0}}}{i\phi_{1}}
$$
より
$$
\begin{align}
C(x)&=\frac{\sin(\phi_{0}+\phi_{1}x)-\sin\phi_{0}}{\phi_{1}} \\
S(x)&=-\frac{\cos(\phi_{0}+\phi_{1}x)-\cos\phi_{0}}{\phi_{1}}
\end{align}
$$
$\phi_{1}=0$なら
$$
\int_0^{x} \exp {i\phi_{0}}\mathrm{d}\sigma=xe^{i\phi_{0}}
$$
より
$$
\begin{align}
C(x)&=x\cos\phi_{0} \\
S(x)&=x\sin\phi_{0}
\end{align}
$$
となります。これで全ての場合の一般化Fresnel積分が求まりました。

ZMでの実装を最後に示します。

/* generalized Fresnel integral */
bool zFresnelIntgGen(double x, double f0, double f1, double f2, double *s, double *c)
{
  double s1, c1, s2, c2;
  bool ret1 = true, ret2 = true;

  if( zIsTiny( f2 ) ){
    zSinCos( f0, &s1, &c1 );
    if( zIsTiny( f1 ) ){
      *s = x * s1;
      *c = x * c1;
    } else{
      zSinCos( f0 + f1 * x, &s2, &c2 );
      *s =-( c2 - c1 ) / f1;
      *c = ( s2 - s1 ) / f1;
    }
  } else{
    ret1 = zFresnelIntgScale( x+0.5*f1/f2, f2, &s1, &c1 );
    ret2 = zFresnelIntgScale(   0.5*f1/f2, f2, &s2, &c2 );
    s1 -= s2;
    c1 -= c2;
    zSinCos( f0 - 0.25 * f1*f1 / f2, &s2, &c2 );
    *c = c1*c2 - s1*s2;
    *s = c1*s2 + s1*c2;
  }
  return ret1 && ret2;
}

クロソイド曲線

せっかくなのでクロソイド曲線を描いてみました。(C(x),S(x))の軌跡をプロットしたものです。
clothoid.png
$\phi_{0}$、$\phi_{1}$、$\phi_{2}$を色々変えることで中心位置、傾き、長さが変わります。

1
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
1