月別アーカイブ: 2016年5月

プログラミングと数値計算, 数学

特殊関数

3件のコメント

以下に示すのは

・ラゲール陪多項式(Associated Laguerre Polynomials):\(L_n^{(\alpha)}(x)\)
・ルジャンドル陪多項式(Associated Legendre Polynomials):\(P_l^m(x)\)
・ヤコビ多項式(Jacobi Polynomials)\(P_n^{(\alpha,\beta)}(x)\)
・エルミート多項式(Hermite Polynomials):\(H_n(x)\)
とその1階微分、2階微分

のfortran90用のプログラムです。
複素平面で有効です。

このプログラムを使用して生じた問題に対する責任は一切取りません。
ミスがあるかもしれないことを念頭に置いて使ってください。
また、使えない場合がありましたらご連絡いただけると幸いです。

※2017/10/27
整数次数の、実数、複素引数エルミート多項式を追加しました。

特殊関数のモジュール”specialfunction”

プログラム例

呼び出す際は総称名
関数：Laguerre(n,a,x), legendre(l,m,x), jacobi(n,a,b,x)
1階微分：Laguerre1(n,a,x), legendre1(l,m,x), jacobi1(n,a,b,x)
2階微分：Laguerre2(n,a,x), legendre2(l,m,x), jacobi2(n,a,b,x)
を用いてください。

ラゲール陪多項式\(L_n^{(a)}(x)\)の引数の型(n,a,x)の組み合わせは、
(n,a,x) : 整数型、　整数型、　　倍精度型　:　出力は倍精度型
(n,a,x) : 整数型、倍精度型、　　倍精度型　:　出力は倍精度型
(n,a,x) : 整数型、倍精度型、複素倍精度型　:　出力は複素倍精度型

ルジャンドル陪多項式\(P_l^m(x)\)の引数の型(l,m,x)の組み合わせは、
(l,m,x) : 整数型、整数型、　　倍精度型　:　出力は倍精度型
(l,m,x) : 整数型、整数型、複素倍精度型　:　出力は複素倍精度型

ヤコビ多項式\(P^{(a,b)}_n(x)\)の引数の型(n,a,b,x)の組み合わせは、
(n,a,b,x) : 整数型、　　倍精度型、　　倍精度型、複素倍精度型　:　出力は複素倍精度型
(n,a,b,x) : 整数型、複素倍精度型、複素倍精度型、複素倍精度型　:　出力は複素倍精度型

というものが有効になっています。

実際に全部使っているプログラムはこちら。

program main
  use specialfunction
  implicit none
 
  write(6,*)laguerre(4,2,4.78d0)
  write(6,*)laguerre(3,2.5d0,4.78d0)
  write(6,*)laguerre(8,dcmplx(1.6d0,3.01d0),dcmplx(5.64d0,49.1d0))
 
  write(6,*)legendre(8,5,0.28d0)
  write(6,*)legendre(6,2,dcmplx(1.8d0,2.7d0))

  write(6,*)jacobi(3,2.15d0,1.85d0,dcmplx(4d0,8.1d0))
  write(6,*)jacobi(3,dcmplx(1.2d0,0.5d0),dcmplx(5d0,0.3d0),dcmplx(4d0,8.1d0))

  stop
end program main

引数は適当に入れています。
spfunc.f90にmoduleを入れ、
temp.f90に上の12行のコードを書きコンパイル→実行すると、

$ ifort spfunc.f90 temp.f90
$ ./a.out    
    3.2997056066666568    
  -8.4458666666671553E-002
 (  454556598.95966828     ,  376437598.22472727     )
   9378.9803534037801    
 ( -519656.04147187521     ,  135464.19972749997     )
 ( -10936.767937499992     , -2092.2704999999974     )
 ( -18248.193121999979     , -10218.566324666655     )
$

と言う結果を得ます。
wolfram alphaと確かめてみましょう。計算で入れた値はそれぞれ、
\(
\begin{align}
& L_4^2(4.78) \\
& L_3^{2.5}(4.78) \\
& L_8^{1.6+i3.01}(5.64+i49.1) \\
& P_8^{5}(0.28) \\
& P_6^{2}(1.8+i2.7) \\
& P_3^{(2.15,1.85)}(4+i8.1) \\
& P_3^{(1.2+i0.5,5+i0.3)}(4+i8.1) \\
\end{align}
\)
です。

wolframの結果と比較したのがこちら。
spfunc_wolfram
下線部が上記プログラムとwolframとで一致した桁です。
まぁまぁいいんじゃないでしょうか。
\(L_8^{1.6+i3.01}(5.64+i49.1)\)の結果が著しく良くないのが気になります･･･使う際は安定な範囲を確かめてから使ってください。

以下、ここで示している特殊関数の定義と関係式です。

ラゲール陪多項式

  • 定義域

    \(\;\;\;\;
    0 \lt x \lt \infty
    \)

  • 微分方程式

    \(\;\;\;\;
    \displaystyle
    \left[x\frac{d^2}{dx^2}+(\alpha+1-x)\frac{d}{dx}+n\right]L_m^{(\alpha)}(x)=0
    \)

  • 具体例

    \(\;\;\;\;
    \begin{align}
    & L^{(0)}_0(x)=1 \\
    & L^{(0)}_1(x)=-x+1 \\
    & L^{(0)}_2(x)=\frac{1}{2}(x^2-4x+2) \\
    & L^{(0)}_3(x)=\frac{1}{6}(-x^3+9x^2-18x+6) \\
    & L^{(1)}_0(x)=1 \\
    & L^{(1)}_1(x)=-x+2 \\
    & L^{(1)}_2(x)=\frac{1}{2}\left(x^2-6x+6\right) \\
    & L^{(1)}_3(x)=\frac{1}{6}\left(-x^3+12x^2-36x+24\right)
    \end{align}
    \)

  • 漸化式

    \(\;\;\;\;
    \begin{align}
    L_n^{(\alpha)}(x)&=\binom{n+\alpha}{n}a_0(x) \\
    &\hspace{1em}a_{m-1}(x)=1-\frac{n-m+1}{m(\alpha+m)}\cdot x\cdot a_m(x) \\
    &\hspace{5em}(m=n,n-1,n-2,\cdots,1,\;\;a_n(x)=1)
    \end{align}
    \)
    [1]
    1階微分
    \(\;\;\;\;
    \begin{align}
    \frac{d}{dx}L_n^{(\alpha)}(x)=-L_{n-1}^{(\alpha+1)}(x)
    \end{align}
    \)
    2階微分
    \(\;\;\;\;
    \begin{align}
    \frac{d^2}{dx^2}L_n^{(\alpha)}(x)=L_{n-2}^{(\alpha+2)}(x)
    \end{align}
    \)

  • 直交性

    \(\;\;\;\;
    \displaystyle
    \int_{0}^{\infty} e^{-x}x^{\alpha}L^{(\alpha)}_n(x)L^{(\alpha)}_{n’}(x) dx =\frac{\Gamma(\alpha+n+1)}{n!}\delta_{n,n’},\;\;(\alpha>-1)
    \)
    [2]

ルジャンドル陪多項式

  • 定義域

    \(\;\;\;\;
    1 \lt x \lt -1
    \)

  • 微分方程式

    \(\;\;\;\;
    \displaystyle (1-x^2)\frac{d^2}{dx^2}P_l^m(x)-2x\frac{d}{dx}P_l^m(x)+\left[l(l+1)-\frac{m^2}{1-x^2}P_l^m(x)\right]=0
    \)

  • 具体例

    \(\;\;\;\;
    \begin{align}
    & P_0^0(x)=1 \\
    & P_{1}^{-1}(x)=\frac{1}{2}(1-x^2)^{1/2} \\
    & P_{1}^{0}(x)=x \\
    & P_{1}^{1}(x)=-(1-x^2)^{1/2} \\
    & P_{2}^{-2}(x)=\frac{1}{8}(1-x^2) \\
    & P_{2}^{-1}(x)=\frac{1}{2}x(1-x^2)^{1/2} \\
    & P_{2}^{0}(x)=\frac{1}{2}(3x^2-1) \\
    & P_{2}^{1}(x)=-3(1-x^2)^{1/2} \\
    & P_{2}^{2}(x)=3(1-x^2)
    \end{align}
    \)

  • 漸化式

    用いている式
    \(\;\;\;\;
    \begin{align}
    & P_{|m|}^{|m|}(x)=(-1)^{|m|}(2|m|-1)!!(1-x^2)^{\frac{|m|}{2}} \\
    & P_{l}^{-m}(x)=(-1)^m\frac{(l-m)!}{(l+m)!}P_l^m(x) \\
    & P_{|m|}^{|m|+1}(x)=(2|m|+1)xP_{|m|}^{|m|}(x) \\
    & P_{|m|}^{|m|+q}(x)=\left(\frac{2|m|-1}{q}+2\right)xP_{|m|}^{|m|+q-1}(x)-\left(\frac{2|m|-1}{q}+1\right)P_{|m|}^{|m|+q-2}(x)
    \end{align}
    \)[3]
    1階微分
    \(\;\;\;\;
    \begin{align}
    m=\pm 1&\; \\
    &:\hspace{1em} \frac{d}{dx}P_l^m(x)=\frac{lx}{(x^2-1)}P_l^m(x)-\frac{l+m}{(x^2-1)} P_{l-1}^m(x) \\
    m\ne\pm 1\; \\
    &:\hspace{1em} \frac{d}{dx}P_l^m(x)=c_{2}P^{m+2}_{l-1}(x)+c_0 P^{m}_{l-1}(x)+c_{-2} P^{m-2}_{l-1}(x) \\
    &\hspace{1em} c_{2}=\frac{1}{4(m+1)}\\
    &\hspace{1em} c_{0}=\frac{l+m}{2}\left(1+\frac{l}{1-m^2}\right)\\
    &\hspace{1em} c_{-2}=-\frac{(l+m)(l+m-1)(l+m-2)(l-m+1)}{4(m-1)}
    \end{align}
    \)
         ※\(m=\pm 1\)の時、\(x=\pm 1\)で発散します。
    2階微分
    \(\;\;\;\;
    \begin{align}
    m=\pm 1,\pm 3&\; \\
    &:\hspace{1em} \frac{d^2}{dx^2}P_l^m(x)=\frac{1}{1-x^2}\left[2x\frac{dP_l^m(x)}{dx}+\frac{(l+1)(l+m)}{2l+1}\frac{dP_{l-1}^m(x)}{dx}-\frac{l(l-m+1)}{2l+1}\frac{dP_{l+1}^m(x)}{dx}\right] \\
    m\ne\pm 1,\pm 3&\; \\
    &:\hspace{1em} \frac{d^2}{dx^2}P_l^m(x)=c_{2}\frac{dP^{m+2}_{l-1}(x)}{dx}+c_0 \frac{dP^{m}_{l-1}(x)}{dx}+c_{-2} \frac{dP^{m-2}_{l-1}(x)}{dx}
    \end{align}
    \)
         ※\(m=\pm 1,\pm 3\)の時、\(x=\pm 1\)で発散します。

  • 直交性

    \(l\)に対する直交性
    \(\;\;\;\;
    \displaystyle
    \int_{-1}^{1} P_{m}^{l}(x)P_{m}^{l’}(x) dx =\frac{2}{2l+1}\cdot \frac{(l+m)!}{(l-m)!}\delta_{l,l’}
    \)
    [2]

ヤコビ多項式

  • 定義域

    \(\;\;\;\;
    -1 \lt x \lt 1
    \)

  • 微分方程式

    \(\;\;\;\;
    \displaystyle
    \left[(1-x^2)\frac{d^2}{dx^2}+(\beta-\alpha-(\alpha+\beta+2)x)\frac{d}{dx}+n(n+\alpha+\beta+1)\right]P_n^{(\alpha,\beta)}(x)=0
    \)

  • 漸化式

    \(\;\;\;\;
    \begin{align}
    P_n^{(\alpha,\beta)}(x)&=\binom{n+\alpha}{n}a_0(x) \\
    &\hspace{1em}a_{m-1}(x)=1-\frac{(n-m+1)(\alpha+\beta+n+m)}{2m(\alpha+m)}(1-x)\cdot a_m(x)\\
    &\hspace{5em}(m=n,n-1,n-2,\cdots,1,\;\;a_n(x)=1)
    \end{align}
    \)
    [1]
    1階微分
    \(\;\;\;\;
    \begin{align}
    \frac{d}{dx}P_n^{(\alpha,\beta)}(x)=\frac{\Gamma(\alpha+\beta+n+2)}{2\Gamma(\alpha+\beta+n+1)}P_{n-1}^{(\alpha+1,\beta+1)}(x)
    \end{align}
    \)
    2階微分
    \(\;\;\;\;
    \begin{align}
    \frac{d^2}{dx^2}P_n^{(\alpha,\beta)}(x)=\frac{\Gamma(\alpha+\beta+n+3)}{2^2\Gamma(\alpha+\beta+n+1)}P_{n-2}^{(\alpha+2,\beta+2)}(x)
    \end{align}
    \)

  • 直交性

    \(\;\;\;\;
    \begin{align}
    \displaystyle
    &\int_{-1}^{1} (1-x)^{\alpha}(1+x)^{\beta}P_{n}^{(\alpha,\beta)}(x)P_{n’}^{(\alpha,\beta)}(x)dx \\
    &\hspace{4em}=\frac{2^{\alpha+\beta+1}}{2n+\alpha+\beta+1}
    \frac{\Gamma(n+\alpha+1)\Gamma(n+\beta+1)}{n!\Gamma(n+\alpha+\beta+1)}\delta_{nm}\\
    &(\alpha>-1,\ \beta>-1)
    \end{align}
    \)    [2]

参考文献

[1] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 22.18, Page 789)
[2] Abramowitz & Stegun (1965), ISBN 978-0486612720 (Numerical Method, 22.2, Page 774,775)
[3] Associated Legendre polynomials -wikipedia

数学

平方根を手計算で簡単に求める方法

2件のコメント

平方根の値を簡単に得る方法です。

もしかしたら、某福〇漫画に出てくるように、\(\sqrt{2}\)を計算しないと命がヤバいみたいな状況に陥るかもしれません。そんな状況を乗り切りましょう。

開平法による筆算とかニュートン法を手でやる方法とかありますが、
覚えるのが大変ですし、計算過程が大変でミスが起こります。

使うべきは
連分数
です。
連分数を使うと平方根の計算が1回の割り算に変わります(確かめのためには2回の割り算が必要)。

早速本題に入りましょう。
\(\;\;\;\;\sqrt{2}\)
を求めてみます。

まず式変形をします。
\(
\;\;\;\;\sqrt{2}=1+(\sqrt{2}-1)
\)
として、整数部分\(1\)と小数部分\((\sqrt{2}-1)\)に分けます。
そして、小数部分に対し、
\(
\;\;\;\;\displaystyle \frac{\sqrt{2}+1}{\sqrt{2}+1}
\)
を掛けますと、
\(
\;\;\;\;\displaystyle \sqrt{2}=1+\frac{1}{1+\sqrt{2}}
\)
になります。ここで、整数部分を左辺に移行しまして、
\(
\;\;\;\;\displaystyle \sqrt{2}-1=\frac{1}{1+\sqrt{2}}
\)
とします。ここで、
\(
\;\;\;\;x=\sqrt{2}-1
\)
と置いて、\(\sqrt{2}=x+1\)を代入すると
\(\;\;\;\;
\begin{align}
x&=\frac{1}{2+x} \\
x&=(2+x)^{-1} \;\;\cdots (1)
\end{align}
\)
なる関係式が得られます。
\(\sqrt{2}=x+1\)ですから、\(x\)が求まればいいわけです。

ここまで来れば後は簡単です。漸化式と見立てて(1)を逐次的に解けばいいのです。
まず、(1)式に\(x(=\sqrt{2}-1)\)の近似値として0という値を入れます。
この値は\(x\)の近似値であれば何でも構いませんが、後々のために整数か分数にしたほうがいいです。

(1)の右辺に代入します。{　}の中身は、そこで打ち切った時に割り算して求めた場合の結果で,アンダーラインは\(x=\sqrt{2}-1=0.414213562372\cdots\)との一致具合を表します。

[1回目]
\(
\displaystyle (2+0)^{-1}=\frac{1}{2}\;\;\;\{\cdots x=0.5\}
\)
[2回目]
\(
\displaystyle (2+\frac{1}{2})^{-1}=\frac{2}{5}\;\;\;\{: x=0.\underline{4}\}
\)
[3回目]
\(
\displaystyle (2+\frac{2}{5})^{-1}=\frac{5}{12}\;\;\;\{: x=0.\underline{41}66\cdots\}
\)
[4回目]
\(
\displaystyle (2+\frac{5}{12})^{-1}=\frac{12}{29}\;\;\;\{: x=0.\underline{41}379\cdots\}
\)
[5回目]
\(
\displaystyle (2+\frac{12}{29})^{-1}=\frac{29}{70}\;\;\;\{: x=0.\underline{4142}857\cdots\}
\)
[6回目]
\(
\displaystyle (2+\frac{29}{70})^{-1}=\frac{70}{169}\;\;\;\{: x=0.\underline{4142}011\cdots\}
\)
[7回目]
\(
\displaystyle (2+\frac{70}{169})^{-1}=\frac{169}{408}\;\;\;\{: x=0.\underline{41421}568\cdots\}
\)
[8回目]
\(
\displaystyle (2+\frac{169}{408})^{-1}=\frac{408}{985}\;\;\;\{: x=0.\underline{414213}131\cdots\}
\)
[9回目]
\(
\displaystyle (2+\frac{408}{985})^{-1}=\frac{985}{2378}\;\;\;\{: x=0.\underline{414213}624\cdots\}
\)
[10回目]
\(
\displaystyle (2+\frac{985}{2378})^{-1}=\frac{2378}{5741}\;\;\;\{: x=0.\underline{4142135}551\cdots\}
\)
[11回目]

\(
\displaystyle (2+\frac{2378}{5741})^{-1}=\frac{5741}{13860}\;\;\;\{: x=0.\underline{41421356}421\cdots\}
\)
…

ただ2倍して足して、分子分母をひっくり返せば計算できるので、何も難しいことはありません。
ここまでの計算は5分もかからないでしょう。

しかも、途中で計算ミスをしても構わないのです。
なぜなら、近似値からさらに良い値にどんどん更新していくからです。計算回数を多くしないといけないだけです。

収束性を確かめたいのなら、もう一回増やしたものを計算する必要があります。その差を見ると良いでしょう。

一回だけ割り算をしなければなりませんね。これはどうしてもなくせません。
ただし、\(\sqrt{2}\)が1.41421356と分かっていれば、筆算の計算は途中までは掛け算と引き算で行けます。そう考えるとまぁ多少は簡単なんですかね・・・？

参考文献

連分数の応用 -講義ノート、筑波大