Subsections

• 1 EXTERNAL文
• 2 問題：Newton-Raphson法を簡単な問題に適用する
• 3 問題：Mollwide画法で天球図を描く

1 Newton-Raphson法

Newton-Raphson法は、関数$f(x)$が単調連続で変曲点がなく、且つ導関数$f'(x)$が 数値的に求められる時に利用できる数値解法で、収束が速く、広く用いられている。 いま、曲線上の点$(x,f(x))$での接線と$x$-軸との交点を$(x_{n+1}, 0)$とすると、

$$f'(x_n)=-{{f(x_n)}\over{x_{n+1}-x_n}}$$ (9)

だから、 解の適当な推定値から始めて、反復公式

$$x_{n+1}=x_n-{{f(x_n)}\over{f'(x_n)}} \qquad (n=0, 1, 2, \cdots) .$$ (10)

を繰り返すと、 漸近的に解に達する。

Figure 6.1: Newton-Raphson法の概念図。

収束の判定条件としては、ある$\varepsilon$に対して、

1. 二つの近似解間の差 $\vert x_{n+1}-x_n\vert<\varepsilon$
2. 二つの近似解間の相対差 $\vert{{x_{n+1}-x_n}\over{x_n}}\vert<\varepsilon$
3. 二つの近似解関数値間の差 $\vert f(x_{n+1})-f(x_n)\vert<\varepsilon$

等が考えられるが、1と2の条件では、$f'(x_n)$の値が大きいと、真の解に未だ 充分近くないのに満たされてしまう場合があるので、3の条件を判定条件として用いる 方が誤った解を出してしまう危険性がすくない。

以下の外部関数rtnewtは、関数$f(x)$をfuncdとして、$f(x)=0$となる解を Newton-Raphson法で解くプログラム例である 。収束判定条件は上記1を使っている。

      function rtnewt(funcd,x1,x2,xacc)
      Implicit NONE
      Integer JMAX
      Real rtnewt,x1,x2,xacc
      external funcd
      Parameter (JMAX=20)       !Set to maximum number of iterations.
!	Using the Newton-Raphson method, find the root of a function known to
!	lie in the interval [x1,x2]. The root 'rtnewt' will be refined until 
!	its accuracy is known within +-'xacc'. 'funcd' is a user-supplied 
!	subroutine that returns both the function value and the first 
!	derivative of the function at the point x.
      Integer j
      Real df,dx,f
      rtnewt=.5*(x1+x2)         !Initial guess.
      do j=1,JMAX
        call funcd(rtnewt,f,df)
        dx=f/df
        rtnewt=rtnewt-dx
        if((x1-rtnewt)*(rtnewt-x2).lt.0.)
     *	   pause 'rtnewt jumped out of brackets'
        if(abs(dx).lt.xacc) return     !Convergence.
      end do
      pause 'rtnewt exceeded maximum iterations'
      end

1 EXTERNAL文

上記サンプルプログラムで、external funcdという行がある。これは、 EXTERNAL文と呼ばれる宣言文である。FUNCTION文の引数は、通常はその関数の独立変 数であるが、今の場合、funcdは、解きたい関数名である。そして、 この関数は、別にサブルーチンで与える。そこで、 FUNCTION文の引数が、外部手続きサブルーチンの名である事を宣言しておくのである。

2 問題：Newton-Raphson法を簡単な問題に適用する

Newton-Raphson法で、

$$f(x)=2x+\sin{x}-2\pi=0$$ (11)

を解け。

3 問題：Mollwide画法で天球図を描く

CfA カタログを使って、銀河の天球上の 赤経・赤緯座標、銀経・銀緯座標、の二次元マップを、 Mollwide画法で描いて見よう。 Mollwide図法では経度・緯度$\lambda$-$\phi$を次の様に$x$-$y$座標に変換する。

$$x=(\sqrt{8}/\pi) (\lambda-\lambda_0)\cos\theta,$$ (12)

$$y=\sqrt{2}\sin\theta.$$ (13)

但し、$\lambda_0$は地図上の中心経度であり、$\theta$は

$$2\theta+\sin 2\theta = \pi\sin\phi$$ (14)

で定義する 。 （経度$\lambda$はファイルGALAXY.DATから読み取った値を直接使用するのではない。 適切な変換をしてから使用する。 図6.2、6.3の$x$-軸を見て、良く考えよう。） 与えられた緯度$\phi$から上の式に基づいて$\theta$を求めるには、 先ず適当な推定値$\theta$を用意して、 Newton-Raphson法により、上の式を満たすように逐次修正して行けば良い。 修正値$\Delta\theta$は、(6.2)式を使って、

$$\Delta \theta=-(2\theta+\sin2\theta-\pi\sin\phi)/2(1+\cos2\theta)$$ (15)

で求め、これを$\theta$に加えたものを改めて$\theta$として、 修正値$\Delta\theta$が予め定めておいた限界値以下になるまで、 上のプロセスを繰り返す。

Figure 6.2: 銀河の天球上の分布。赤経、赤緯座標で図示。

赤経($\alpha$)・赤緯($\delta$)と銀経($l$)・銀緯($b$)の変換は、 銀河の北極の赤経・赤緯を$\alpha_0$，$\delta_0$とし、

$$\Omega=90^\circ + \alpha_0$$ (16)

$$i = 90^\circ - \delta_0$$ (17)

$$\alpha' \equiv \alpha - \Omega$$ (18)

として、

$$\cos b \cos l = \cos\delta\cos\alpha'$$ (19)

$$\sin b = \cos i \sin\delta - \sin i \cos\delta\sin\alpha'$$ (20)

$$\cos b \sin l = \sin i \sin\delta + \cos i \cos\delta\sin\alpha'$$ (21)

で与えられる。但し、 $(\alpha_0, \delta_0) = (12^{\rm h}51^{\rm m}.4, +27^{\circ}08'.0)$。

Figure 6.3: 銀河の天球上の分布。銀経、銀緯座標で図示。