数値天文学入門
ー天文学で用いる数値技法ー
福島登志夫
東京大学、総合研究大学院大学
2006
目次









1. 計算の基本
2. 特殊関数
3. 楕円関数
4. 補間と外挿
5. 関数の近似
6. 連立1次方程式
7. 非線形方程式
8. 最適化問題
9. 最小2乗法







10. 数値微分
11. 数値積分
12. 常微分方程式の数
値解法 I. 1段法
13. 常微分方程式の数
値解法 II. 多段法
14. 常微分方程式の数
値解法 III. その他
15. 計算のテクニック
16. 計算の高信頼性化
はじめに

ある天文学者の日常生活
メールの読み書き、情報検索、文書作成
 機器制御、データ管理、データ解析
 シミュレーション、プレゼン作成、論文執筆


数値天文学とは

(天文学で行う)データ解析やシミュレーショ
ンに必要な数値計算技法
はじめに(2)

本講義で扱わない分野


偏微分方程式、積分方程式、微分・代数方
程式、統計、乱数、ソート…
参考書
名著:Ralston and Rabinowitz (邦訳あり)
 名著者:Kernighan(邦訳書多数)
 網羅的だが玉石混交: Numerical Recipes


個人用ならフリーでダウンロード可
1. 計算の基本







基本演算のコツ
浮動小数点演算とは
数学ソフトウェアと数学ライブラリ
定義式による計算
テイラー展開、パデ近似、連分数
基本変数から出発
中間結果の再利用、漸化式
基本演算のコツ


基礎知識:浮動小数点数の表現法
常に倍精度で計算するのが良い
ただし、メモリー不足にならない限り
 もし使える環境にあるなら拡張精度で


IEEE754浮動小数点演算:拡張精度計算
単精度は、倍精度より少し遅い
 4倍精度は、倍精度より数~数十倍遅い
 複素数は、実数のペアより遅い

浮動小数点数とは



整数(integer)、実数(real)、複素数(complex)
固定(fixed)小数点数:絶対精度重視
浮動(floating)小数点数:相対精度重視


有限桁で効率的に数を近似表現
浮動小数点表現
x  s  mb
e
符号(sign) s=1、基数(base) b=2 or 16
 仮数(mantissa) m、指数(exponent) e

浮動小数点演算の精度
精度 Fortran
C
仮数+指数 マシン・エプシロン
(bit数)
24+8
単
Real*4
float
倍
Real*8
double 53+11
拡張 Real*10 4倍

64+16
Real*16 long
113+15
double
2-24 ~ 5.96x10-8
2-53 ~ 1.11x10-16
2-64 ~ 5.42x10-20
2-113 ~ 9.63x10-35
ケチ表現による仮数オマケ1ビット & 符号1ビット
スケーリングの重要性


浮動小数点表現を最大限に活用
教訓:絶対値が大きい数値を扱わない


要点:1程度の量に変換して計算


悪例:惑星の軌道半径をmで記述
オーバーフローの回避(特に大型行列など)
手法:無次元化=代表(nominal)値で除算
=x/x0, u=v/v0など
 丸め誤差の防止:代表値は2の累乗が良い

数学定数の計算
なるべく組み込まれた数学関数を使う
e : exp(1.0),  : 4.0*atan 1.0  ,
2 : sqrt  2.0 

頻繁に使う定数は、 1
n 
n
, 2 , n!, n!!,  
1回だけ計算して保存 n
 m
 超越定数の近似値
  3.1415926535897932384626433832795

e  2.7182818284590452353602874713527
数学ソフトウェア

問:調査せよ(参考:山本他1999)
Mathematica, Maple: 純粋数学志向、高価
 MATLAB: 行列に特化、応用数学志向、高価
 MAXIMA: フリー、Mathematicaに類似
 scilab: フリー、MATLABに類似
 awk, bc, perl: そっけないが使い良い
 GnuMP: フリー、任意多倍長演算
 GNUPLOT: フリー、グラフ作成

数値計算ライブラリ


定番:IMSL, NAG, LAPACK
日本発:NUMPAC



http://netnumpac.fuis.fukui-u.ac.jp/numpac/
手軽だが要注意:Numerical Recipe
無料(ただ)ほど便利なものは無い
NetLIB: なんでもそろうデパート
 GSL(=GNU Sci. Lib.): C, C++
 SLATEC: Fortran, FFTW: FFT専用

定義式による計算


非効率なことが多い
Cnm   Ank Bkm
効率的な例:行列の積
k


例外:シュトラッセン(Strassen)法
非効率な例
桁落ちが起きる場合: x~0での1-cos(x)
 逆行列を用いた連立1次方程式の解
 多項式、連分数、フーリエ級数
 直交多項式、球面調和関数など

テイラー(Taylor)展開

多項式近似の基礎

f  x  
微分が容易なら有効
n 0
 多項式の計算はホーナー法で


f
 n
 x0 
n!
 x  x0 
n
いくつかの例(注:cnはn次シュトゥンプ関数)
cosh x  c   x 
x 

x
x
cos x  
 1    c  x 
2 24
 2n  !

2 n
2
2
4
2
0
n 0

0
x 
2 n
x3 x5
sin x  x
 x 

6 120
n  0  2n  1 !
 xc1  x 2 
sinh x  xc1   x 2 
e x  cosh x  sinh x
ホーナー(Horner)法

多項式の計算の定番

定義式で計算:計算量=N(N+3) /2



P  x    an x n
n 0
計算量=加減算と乗算の数~Flops
入れ子で計算:計算量=2N-1、丸め誤差も低減
P  x   a0  x  a1  x  a2  x  a3    
アルゴリズム:添字の向きに注意(問:なぜ)
P:=a N ;do  k=N-1,0,-1P:=a k +x*P

N
ベクトル多項式:さらに工夫が必要
2次因子法

多項式計算の決定版


単純2次式の入れ子で計算:計算量=1.5N+1
2
P  x    y  c0   b1   y  c1   b2   y  c2  b3    
y   x  a
欠点:準備に多大な計算コスト
理由:多数の高次代数方程式を解く必要
 詳細:Ralston and Rabinowitz 7.2
 頻繁に用いる関数等でない限り不必要

シュトゥンプ(Stumpff)関数


Stumpff (Himmelsmechanik 1-3, 1959-74)
三角関数と双曲線関数の統一的扱い
z 

1
z
z2
cn  z   
 


n !  n  2 !  n  4 !
k  0  n  2k  !

2
 z=x >0
k
1
  zcn  2  z 
n!
c0  x 2   cos x, c1  x 2    sin x  / x, c2  x 2   1  cos x  / x 2
2
2
2
2
2
c

x

cosh
x
,
c

x

sinh
x
/
x
,
c

x

cosh
x

1
/
x










0
1
2
 z=-x <0

一般調和振動解

問:確認せよ
d2 x
2
2


x

0

x

x
c

t

v
tc

t




0 0
0 1
2
dt
シュトゥンプ関数(2)



正規化シュトゥンプ関数 Cn  z   n!cn  z 
Cn  0  1
逆行漸化式 Cn  1 n  n 1 zCn2
実用的な計算法(問:コードを書け)
 1. 初期値:十分大きいN CN  CN 1  1

N=8(|z|<0.01), 10(|z|<0.1), 19(|z|<1)
2. 逆行漸化式(添字が偶・奇で並列計算)
Cn
 3. (必要なら)標準形への変換
cn 
c0  C0 , c1  C1
n!

テイラー展開(2)

いくつかの例(続き)

詳細は超幾何関数表現から
x3 2 x5 17 x7 62 x9 1382 x11 21844 x13 929569 x15
tan x  x  






3 15 315 2835 155925 6081075 638512875
x3 3x5 5 x 7 35 x9 63x11 231x13 143x15
sin x  x  






6 40 112 1152 2816 13312 10240
1
xn
x 2 x3
 log 1  x     x   
2 3
n 1 n

tan 1 x  x
n 0
x 
x3 x5
 x  
2n  1
3 5
 xs1  x 2 
tanh x  xt   x 2 

2 n
 xt  x 2 
 xt1  x 2 
sinh 1 x  xs1   x 2 
tanh 1 x  xt1   x 2 
ガウスの超幾何関数

超幾何(hypergeometric)関数




多くの初等関数の統一的扱い
級数による定義
  a n  b n  x n
F  a, b, c; x    

n 0 
  c n  n !

ポッホハンマー(Pochhammer)記号
  a  n
Fの性質
 a  a  1  a  n  1
a 
 a
n
 
aとbについて対称
 a or bが負整数なら有限和→超幾何多項式

合流型超幾何関数


クンマー(Kummer)
合流型(congruent)超幾何関数


別種の特殊関数の統一的扱い
級数による定義

  a n  x n
F  a, c; x    

n 0 
  c n  n !

aが負整数→合流型超幾何多項式
一般化超幾何関数



ポッホハンマー(Pochhammer)
一般化(generalized)超幾何関数
定義 F a , , a ; b , , b ; x     a1 n  a p n  xn
   b  b n!
p q 1
p 1
q
n 0 
 1 n  q n 
a1, a2, …, apのどれかが負整数→多項式
 問:ホーナー流のpFq計算プログラムを書け


ガウス超幾何関数と合流型超幾何関数
F  a, b, c; x   2 F1  a, b; c; x  F  a, c; x   1 F1  a; c; x 
超幾何関数の応用

初等関数の超幾何関数表現
e x  0 F0  x  ,
1  x 
a
 1 F0  a; x  ,
 3  x2 
sin x  x 0 F1  ;
,
2 2 
 log 1  x   xF 1,1, 2; x  ,
 3 x2 
sinh x  x 0 F1  ;  ,
2 2 
 1  x2 
 1 x2 
cos x  0 F1  ;
 , cosh x  0 F1  ;  ,
2 2 
2 2 
1 1 3 2
1 1 3

1
1
sin x  xF  , , ; x  , sinh x  xF  , , ;  x 2  ,
2 2 2

2 2 2

1 3

tan x  xF  ,1, ;  x 2  ,
2 2

1
1
 1 x 
1 3 2
tanh x  log 
  xF  ,1, ; x 
2
 1 x 
2 2

1
有理式近似

有理式近似の二大手法


パデ(Padé)近似と連分数(continued fraction)
有理式近似の例
ex  1 
x
x

1


Bk x k
x x2 x4
1  


2
12
720
k
!
k 0
x3 x5 97 x 7
x  

12 30 5040
log 1  x  
x x3 7 x5 181x7
1  


2 12 180 7560
問:テイラー展開と誤差・計算速度を比較せよ
 Bn:ベルヌーイ(Bernouli)数、 B2n+1=0 (n>1)

 1
n

 1 1 1 1 5 691 7 3617 43867 174611 854513 236364091
B2 n   , , , , ,
, ,
,
,
,
,
,
330
138
2730
 6 30 42 30 66 2730 6 510 798



パデ近似

(m,n)次パデ近似

分子m次式、分母n次式、近似次数=m+n
m
Rmn  x  
k
p
x
 k
k 0
n
1  q j x j
 f  x   O  x m  n 1 
j 1


m=n-1 or n → m+n=一定の中で最も精確
計算コスト:m+n次テイラー展開と互角

主要部分に偶奇性→より高速(~倍速)
パデ近似(2)

テイラー展開からの構成法




m n
f  x    fi xi
i 0
仮定:mとnは固定
分子多項式 p(x) p  x   m p x k q  x   1  n q x j

j

k
j 1
k 0
分母多項式 q(x)
係数決定条件 p  x   f  x  q  x   O  x 
m  n 1

決定(連立1次)方程式の解→q、決定条件→p
n
f
j 1
m i  j
q j   f m i
 i  1,
, n
k
pk  f k   f k  j q j
j 1
パデ近似(3)

パデ近似の特徴
長所:同次数テイラー展開より誤差係数が小
 欠点:次数の変更→係数の再計算が必要


一例:指数関数の(n,n)次パデ近似

問:テイラー展開と誤差係数・計算速度を比べよ
x
1
x3
2
exp  x  
 
x 12
1
2
x x2
1 
x5
2
12



2
x x
720
1 
2 12
x x 2 x3
1  
x7
2
10
120



2
3
x x
x
100800
1  
2 10 120
連分数展開


連分数
an |
a1
b0  
 b0 
n 1 | b
a2
n
b1 
b2 
定義と略記法
 長所:テイラー展開より収束が良い
 欠点:除算の連続→計算コスト高


連分数級数の例:tan(x)
問:打ち切り次数を
変えて、テイラー展開と
誤差、速度を比較せよ

x
x
tan x 

2

x
x2
1 
1
n 1 2 k  1
x2
3
5
連分数展開(2)

連分数展開の実用計算法:漸化式の活用
発想:単純分数式に
変換後、除算を1回
 結果~パデ近似



ak | An


k 1 | b
Bn
k
n
最初の数項 A0  0, B0  1, A1  a1 , B1  b1
三項漸化式 An  bn An1  an An2

問:導け
Bn  bn Bn1  an Bn2
パデ近似(4)

別の例:tan(x)の連分数展開を整理
x3
x
x
5
15  O x 7
tan x 

O
x





2
2
x
2x
1
1
3
5
3
5
3
x
x
2x
x 
x
9
9 945  O x11
21


O
x





2
4
2
4
3x
x
4x
x
1

1

7 105
9
63
変数変換による加速

発想:変数変換→収束が遅い級数の加速
a  bx
適切な変換を発見するのは職人芸
y
 よくある例:一次有理変換
c  dx
n

例:対数関数
x

x 2 x3
log 1  x   
 x  
2 3
 テイラー展開
n 0 n  1



x>0では1次交代級数→収束が遅い
一次有理変換
 加速後


x
y
2 x

 1 y 
y 2n
log 1  x   log 
  2 y
2次単調増加級数
n 0 2n  1
 1 y 
変数変換による加速(2)

k

k  
1
オイラー(Euler)変換 a 
a


n
k 1     j 
n 0
k 0 2
 j 0  j  
 ゆっくり変化する交代級数の計算に威力


a0 b0  1  k 1  k  1 
an  an 1  bn   an     k 1   
bj 

2 4 k  2 2  j 0  j  
n 0


問:次の無限級数を4桁精度で求めよ


ヒント:桁落ちに注意してbnを書き換えよ
応用:振動積分の計算
an 
 1
 f  x  sin  kx  dx
3
n 1
n
基本変数から出発

三角関数:基本は半角正接関数

t  tan

同じ引数でsinとcosを計算する場合に好都合

2
2t
1 t2
2t
sin  
, cos  
, tan  
2
2
2
1 t
1 t
1 t
双曲線関数:基本は半角正接双曲線関数
x
2t
1 t2
2t
t  tanh sinh x 
, cosh x 
, tanh x 
2
2
2
1 t
1 t
1 t2
例:sinとcosの同時計算

同じ引数のsinとcosの双方を計算

仮定:0<</4
t:=tan  θ*0.5 ; t2:=t*t; d:=1.0/(1.0+t2)
s:=(t*2.0)*d; c:=(1.0-t2)*d
 問:上記アルゴリズムを用いて、任意角度の
sinとcosの同時計算ルーチンsincosを書け
 問:数学ライブラリのsinとcosを、それぞれ呼
ぶ場合と計算誤差・計算速度を比較せよ
基本変数から出発(2)

atan2
y
,
x


逆三角関数:基本は2変数逆正接


sin 1 s  atan2 s, 1  s 2 ,cos1 c  atan2



1  c2 , c , tan 1 t  atan2 t,1
逆双曲線関数と対数関数

tanh-1の引数が小のとき:基本は逆正接双曲線
s
sinh 1 s  tanh 1

1  s2
c2 1
,
c
, cosh 1 c  tanh 1
 x 1 
log x  2 tanh 1 

 x 1 
logの引数が大のとき:基本は対数




sinh 1 s  sign  s  log s  1  s 2 , cosh 1 c  log c  c 2  1 ,
1
 1 t 
tanh 1 t  log 

2
 1 t 
中間結果の再利用

必要な中間結果は、保存して何度でも使う
例1: 同じ除数による除算→逆数による乗算
 例2: ベクトル級数
do  k=0,K  c k :=F  θ k 

f   Ak F k 
k

do  n=1,N  {f n :=A n0 *c0
do  k=1,K f n +=A nk *c k 
特にベクトル多項式

F(xk)=xk の計算にホーナー法を併用
c0:=1.0; c1:=x; do  k=2,K ck :=ck-1*x
2. 特殊関数



超幾何関数による統一表現
チェビシェフ(Chebyshev)多項式
ルジャンドル(Legendre)多項式




ルジャンドル陪(associated)関数
整数次ベッセル(Bessel)関数
エルミート(Hermite)多項式
ラゲル(Laguerre)多項式
超幾何関数表現

直交関数・直交多項式の超幾何関数表現
負整数パラメータ → 多項式
 チェビシェフ多項式


1 1 x 

Tn  x   F  n, n, ;

2 2 

Tchebysheffとも綴る
ルジャンドル多項式
1 x 

Pn  x   F  n, n  1,1;

2 
 ルジャンドル陪多項式

m
n  m !

1 x 

m
2
Pn  x   m
1 x
F  m  n, n  m  1, m  1;

2 m! n  m!
2





超幾何関数表現(2)

合流型 or 一般化超幾何関数表現

整数次ベッセル関数

エルミート多項式


偶数次

奇数次
ラゲル多項式
1 x
Jn  x   
n!  2 
n

 x2 

0 F1  ; n  1;
4


2


1
x
n
H 2 n  x    1  2n  1!! F  n; ; 
2 2


3 x2 
H 2 n1  x    1  2n  1!! xF  n; ; 
2 2

n
Ln  x   F  n;1; x 
k !!  k  k  2 k  4
漸化(recurrence)式

中間結果の再利用の典型

超幾何関数の三項漸化式のいくつか
bFa ,b1,c   b  a  Fa,b,c  aFa 1,b,c
 c  b  Fa,b1,c   a  b 1  x  Fa,b,c   c  a  Fa1,b,c
 c  a  Fa1,b,c  c  2a   a  b  x  Fa,b,c  a 1  x  Fa 1,b,c

合流型超幾何関数の三項漸化式の一つ
aFa1,c   x  2a  c  Fa,c   c  a  Fa1,c
チェビシェフ多項式

x  cos 
変形三角多項式
第1種T、第2種U T  x   cos n ,U  x   sin n
n
n
sin 
自明な性質



1. ゼロ点
 0x
 n
Tn xk
 n
k

2. 極値点

3. 極値の交代性 Tn  yk
  2k  1  
 cos 

 2n 
dTn  n
 k 
 n
yk  0  yk  cos 

dx
n


 
 n
   1
k
 k  1,..., n 
 k  0,1,..., n 
チェビシェフ多項式(2)
1
T1
0.5
T2
0
T3
-0.5
-1
-1
-0.5
0
0.5
T4
T5
1
チェビシェフ多項式(3)

計算法:漸化式が最良


最初の数項
T0  1,U0  0, T1  x,U1  1
二項漸化式(加法定理):誤差が増大
Tn 1  xTn  1  x U n ,U n 1  xU n  Tn
2

三項漸化式(積和公式):高速、独立計算可
Tn1  2xTn  Tn1,Un1  2xUn  Un1
チェビシェフ多項式(4)

チェビシェフ多項式の変形: Wn, Vn
y~0のときTnの漸化式計算は桁落ちしやすい
 UnよりVnのほうが実用的
Wn  1  Tn ,Vn  yUn  sin n y  sin 

三項漸化式に基づくアルゴリズム(問:導け)
t:=tan  θ*0.5 ; y:=t*2.0/ 1.0+t*t  ; w:=t*y; w2:=w*2.0; z:=r2+2.0;

W0 :=0.0; T0 :=1.0; V0 :=0.0; W1:=w; T1:=1.0-w; V1:=y;
do  n=1,N Wn+1:=z*Wn +w2-Wn-1; Tn+1:=1.0-Wn ; Vn+1:=z*Vn -Vn-1
チェビシェフ多項式(5)

具体的表現(低次の場合) T0  x   1, T1  x   x
T2  x   2 x 2  1
T3  x   4 x3  3x
T4  x   8 x 4  8 x 2  1
T5  x   16 x 5  20 x 3  5 x
T6  x   32 x 6  48 x 4  18 x 2  1
T7  x   64 x 7  112 x5  56 x3  7 x
T8  x   128 x8  256 x 6  160 x 4  32 x 2  1
T9  x   256 x9  576 x 7  432 x5  120 x 3  9 x
チェビシェフ多項式(6)


逆表現=単項式をチェビシェフ多項式で
n / 2 n
 
1
n
x  n1    Tn2k  x 
2 k 0  k 
例 1  T0
16 x 5  T5  5T3  10T1
x  T1
32 x 6  T6  6T4  15T2  20T0
2 x 2  T2  T0
64 x 7  T7  7T5  21T3  35T1
4 x 3  T3  3T1
128 x8  T8  8T6  28T4  56T2  70T0
8 x 4  T4  4T2  6T0 256 x 9  T9  9T7  36T5  84T3  126T1
多重チェビシェフ多項式

チェビシェフ多項式の多次元版

係数(=添字)が整数ベクトル、角度もベクトル
Tn  x  cos n θ ,Vn  x  sin n θ  xk  xk  cosk

摂動論で多用


好例:章動理論(十数次元、数千項)
効率的計算法→三角関数の利用を最小化
(1次元の)漸化式+加法定理
 問:アルゴリズムを導け

ルジャンドル多項式
n


定義 P  x   1  d   x2  1n
n
 
n
2 n!  dx 
性質
1. 対称性
Pn
 2. 特別な点での値

  x    1
n
Pn  x 
2n  1!!

P2 n  0    1 pn 0 , P2 n 1  0   0 pn0 
 2n !!
n
Pn 1  1, Pn  1   1 N !!  N  N  2 N  4
n
ルジャンドル多項式(2)

級数表現
P2 n  x    1
n
 p x 
n
2 k
nk
k 0
2n 

2 k
P2 n1  x    1 x 1 
 pnk   x 
2k  1 
k 0 
n

n
係数の漸化式(問:定義式より導け)
p0,0  1
pn 1,0
 2n  1 

 pn,0
 2n  2 
n  k  2n  2k  3

pn,k 1 
pn,k
 k  1 2k  1
ルジャンドル多項式(3)

具体的表現(低次の場合) P0  x   1, P1  x   x
P  x    3x  1 / 2
P  x    5 x  3x  / 2
P  x    35 x  30 x  3 / 8
P  x    63x  70 x  15 x  / 8
P  x    231x  315 x  105 x  5  /16
P  x    429 x  693x  315 x  35 x  /16
P  x    6435 x  12012 x  6930 x  1260 x  35  /128
2
2
3
3
4
2
5
3
4
5
6
4
2
7
5
3
6
7
8
8
6
4
2
ルジャンドル多項式(4)
1
P1
0.5
P2
P3
P4
0
-0.5
-1
-1
P5
-0.5
0
0.5
1
ルジャンドル多項式(5)

同一引数に対するPnの一斉計算法
三項漸化式  n 1 Pn1   2n 1 xPn  nPn1
 級数表現+ホーナーの方法より誤差が小


アルゴリズム
P0:=1.0;P:=x;do
 n=1,NPn+1:=An *x*Pn -Bn *Pn-1
1
 コツ:数定数は予め計算
An  1  Bn , Bn  1  I n1
 整数の逆数Inは重宝
1
In 
→ 常備すべし
n
ルジャンドル多項式(6)

Pnの根(対称性より正値のみ表示)
数値積分で有用
 注:P6以降は17桁正しい近似値

P3 : 3/ 5, P4 :
15 

120 / 35, P5 :
35 

280 / 63
P6 : 0.23861918608319691, 0.66120938646626451, 0.93246951420315203
P7 : 0.40584515137739717, 0.74153118559939444, 0.94910791234275852
P8 : 0.18343464249564980, 0.52553240991632899,
0.79666647741362674, 0.96028985649753623
ルジャンドル多項式(7)

Pnの根(続) P : 0.32425342340380893,
9
0.61337143270059040,
0.83603110732663579, 0.96816023950762609
P10 : 0.14887433898163121, 0.43339539412924719, 0.67940956829902441,
0.86506336668898451, 0.97390652851717172
P11 : 0.26954315595234497, 0.51909612920681182, 0.73015200557404932,
0.88706259976809530, 0.97822865814605699
P12 : 0.12523340851146392, 0.36783149899818019, 0.58731795428661745,
0.76990267419430469, 0.90411725637047486, 0.98156063424671925
P13 : 0.23045831595513479, 0.44849275103644685, 0.64234933944034022,
0.80157809073330991, 0.91759839922297797, 0.98418305471858815
ルジャンドル多項式(8)

1階微分



dPn
1  dPn 
Qn 



緯度角
dx cos   d 
x  sin 
非球対称ポテンシャルの偏微分に必要
級数表現
Q2 n  x    1
n 1
Q2 n1  x    1
2 x kpnk   x
n
k 1
n

2 k 1
  2n  2k  1 p   x 
n
k 0
2 k
nk
ルジャンドル多項式(9)

1階微分の漸化式
 n 1 Qn1   2n 1 Pn  xQn   nQn1



問:Pnの漸化式を微分して導け
最初の2項 Q0  0, Q1  1
アルゴリズム
Q0 :=0.0; Q1:=1.0
do  n=1,N  Q n+1:=A n *  Pn +x*Q n  -Bn *Q n-1
ルジャンドル多項式(10)

具体的表現(低次の場合) Q0  x  0, Q1  x   1
Q2  x   3 x
Q3  x   15 x 2  3 / 2
Q4  x    35 x 3  15 x  / 2

性質
Q2n  0  0
Q5  x    315 x 4  210 x 2  15  / 8
Q6  x    693 x 5  630 x 3  105 x  / 8
Q7  x    3003x 6  3465 x 4  945 x 2  35  /16
Q8  x    6435 x 7  9009 x 5  3465 x 3  315 x  /16
ルジャンドル多項式(11)

Qnの根(端点を除く正値点のみ表示)
数値積分・常微分方程式の数値解法で有用
 注:Q7以降は17桁正しい近似値

Q3 : 1/ 5, Q4 : 3/ 7, Q5 :
7 

28 / 21, Q6 :
15 

60 / 33
Q7 : 0.20929921790247887, 0.59170018143314230, 0.87174014850960662
Q8 : 0.36311746382617816, 0.67718627951073775, 0.89975799541146016
Q9 : 0.16527895766638702, 0.47792494981044450,
0.73877386510550508, 0.91953390816645881
ルジャンドル多項式(12)

Qnの
根(続)
Q10 : 0.29575813558693939, 0.56523532699620501,
0.784483473663144419, 0.93400143040805913
Q11 : 0.13655293285492755, 0.39953094096534893, 0.63287615303186068,
0.81927932164400668, 0.94489927222288222
Q12 : 0.24928693010623999, 0.48290982109133620, 0.68618846908175743,
0.84634756465187232, 0.95330984664216391
Q13 : 0.11633186888370387, 0.34272401334271285, 0.55063940292864706,
0.72886859909132614, 0.86780105383034725, 0.95993504526726090
Q14 : 0.21535395536379424, 0.42063805471367248, 0.60625320546984571,
0.76351968995181520, 0.88508204422297630, 0.96524592650383857
ルジャンドル陪関数



定義
m
n
P

 x 
1  x2

m
2n n !
d 
 
 dx 
nm
 x 1
2
n
Pn0  Pn
y  cos 
最初の数項 P  1, P  x, P  y
漸化式(多数のうち、最も実用的な組合せ)
m
m1
m1
m
Pm1   2m 1 yPm Pm1   2m 1 yPm
0
0
 n  m 1 P
m
n1

0
1
1
1
  2n 1 xP   n  m P
m
n
注目:Pnの漸化式(m=0)を含んでいる
m
n1
ルジャンドル陪関数(2)

アルゴリズム:数定数は予め計算&保存
In 
1
, Cn  2n  1, Anm  Cn I n m 1 , Bnm   n  m  I n m 1
n
P00 :=1.0; P10 :=x; P11:=y;
do  m=1,M-1 {
m+1
m+1
P
m
m
m
m+1
:=C m *y*P ; P
do  n=m+1,N-1 {
m
n+1
m
n
m-1
m
:=C m *y*P
m
n-1
P :=A nm *x*P -Bnm *P }
;
ルジャンドル陪関数(3)

具体的表現(低次の場合)
P21  3xy, P22  3 y 2
P31  15 x 2  3 y / 2, P32  15 xy 2 , P33  15 y 3
P41   35 x3  15 x  y / 2, P42  105 x 2  15  y 2 / 2, P43  105 xy 3 , P44  105 y 4
P51   315 x 4  210 x 2  1 y / 8, P52   415 x 3  105 x  y 2 / 2,
P53   945 x 2  105  y 3 / 2, P54  945 xy 4 , P55  945 y 5
P61   693x5  630 x3  105 x  y / 8, P62   3465 x 4  1890 x 2  105  y 2 / 8,
P63   3465 x 3  945 x  y 3 / 2, P64  10395 x 2  945  y 4 / 2,
P65  10395 xy 5 , P66  10395 y 6
ルジャンドル陪関数(4)

緯度角による微分



m
n
dP
R 
d
第2種=Qnmより実用的
m
n
Rn0  yQn
最初の数項 R00  0, R10  y, R11  x
漸化式(問:Pnの漸化式の微分より導け)
Rmm11   2m  1  yRmm  xPmm 
Rmm1   2m  1  yRmm 1  xPmm 1 
 n  m  1 R
m
n 1
  2n  1  xR  yP
m
n
m
n
  n  m R
m
n 1
ルジャンドル陪関数(5)

アルゴリズム
R 00 :=0.0; R 10 :=y; R 11:=-x;
do  m=1,M-1 {
m
m
R m+1
:=C
*
y*R
-x*P
m+1
m 
m
m 
R
m
m+1
:=C m *  y*R
m-1
m
m-1
m
-x*P

do  n=m+1,N-1 {
m
R mn+1:=A nm *  x*R nm +y*Pnm  -Bnm *R n-1

球面調和関数展開

例:天体(地球など)の重力ポテンシャル


a
U  r ,  ,    1    
r  n 2  r 
n

P  sin   Cnm cos m  Snm sin m  

m0

n
m
n
問1:極座標(r,,)による偏微分を求めよ
 問2:チェビシェフ多項式Tn,Vnとルジャンドル
陪多項式Pnm,Rnmを使って、ポテンシャルお
よび偏微分を表現せよ
 問3:上記を計算するプログラムを書け

整数次ベッセル関数

定義


 x
Jn  x   
2
テイラー展開
平方恒等式
 1  x 

 
k
!
n

k
!

 2
k 0
k
n 

2k
 J 0  x    2  J n  x   1
2
2
n 1


三項漸化式
微分漸化式
xJ n1  2nJ n  xJ n1
dJ 0
dJ n n
1
  J1
 J n  J n 1
dx
dx x
2
整数次ベッセル関数(2)
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
0
J0
J1
J2J3
J4
2
4
6
8
10
整数次ベッセル関数(3)

補助関数Knの導入
2

1 x
x 
Jn  x     Kn 

n!  2 
 4 
n

k
z
K n  z   0 F1  n  1; z   
k  0  n  1 k k !

補助関数の三項(逆行)漸化式
z
K n1  K n 
K n 1
n  n  1
Kn  0  1
整数次ベッセル関数(4)



ベッセル関数の計算法
1. 小さい引数(|x|<0.1)→テイラー展開
2. 中程度の引数(|x|<n) でJ0とJ1の関数
ライブラリが利用可能
 2n 


J n1    J n  J n 1
J2以降は三項漸化式で
 x 
3. それ以外

漸化式による逆行計算法
整数次ベッセル関数(5)

漸化式による逆行計算法(|x|>0.1)

1. 初期値:十分大きいN jN  0, jN 1  
N=13 (x<1), 30 (x<10), 150 (x<100)
 マシン・エプシロン
 2n 



jn1    jn  jn1
2. 逆行漸化式
 x 
3. 平方恒等式による定数調整
N
jk
2
2
Jk 
C   j0   2  jn 
C
n 1
整数次ベッセル関数(6)

1階微分(n>0)の計算法
dJ n  n 
引数が小さくないとき
1
   J n    J n1
(|x|>0.1):微分漸化式 dx  x 
2
 引数が小さいとき(|x|<0.1)

微分漸化式は桁落ちが激しい
 補助関数による表現(問:示せ)

dJ n
1
 x

 
dx 2  n  1!  2 
n 1
   x2 
  x2 
 Kn 
  K n 1 

 4 
  4 
円盤調和関数展開

例:円盤(銀河など)上の重力ポテンシャル
  

 n
U   ,     J 0     J n     Cnm cos m  Snm sin m 
  a  n1  a  m0

問1:極座標(,)による偏微分を求めよ
 問2:チェビシェフ多項式Tn, Vnと整数次ベッセ
ル関数Jn, J’nを使って、ポテンシャルおよび偏
微分を表現せよ
 問3:上記を計算するプログラムを書け

エルミート多項式



定義
x
n
H n  x    1 exp 
 2
2
 d  
  x2  
    exp 

  dx  
 2 
n
n / 2
級数表現
 n  n2 k
k
H n  x     1  2k  1!!  x
性質
k 0
 2k 
n
1. 対称性
H n   x    1 H n  x 
 2. 特別な点での値

H 2 n  0    1  2n  1!!, H 2 n 1  0   0
n
エルミート多項式(2)
1.2
H0
1
0.8
H5
0.6
H1
0.4
H4
0.2
0
-0.2
-0.4
-0.6 H2 H3
0 0.5 1 1.5 2
  x2 
1
exp 
 Hn  x
n
 2 
2.5
3
3.5
4
エルミート多項式(3)


漸化式 Hn1  xHn  nHn1 Hn1   n 1 Hn
具体的表現(低次の場合)
H0  x  1
H1  x   x
H 2  x   x2 1
H 3  x   x3  3x
H 4  x   x4  6x2  3
H 5  x   x5  10 x3  15 x
ラゲル多項式

定義

級数表現

漸化式

特別な点での値 Ln  0  1, Ln  0  n
n
e  d  x n
Ln  x      e x 
n!  dx 
x
k
n


x
k
Ln  x     1  
k 0
k  k!
n
 n 1 Ln1   2n 1 x  Ln  nLn1
ラゲル多項式(2)
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
0
 x 
exp   Ln  x 
 2 
L0
L2
L5
L4
L3
L1
2
4
6
8
10
ラゲル多項式(3)

具体的表現
(低次の場合)
L0  x   1
L1  x   1  x
x2
L2  x   1  2 x 
2
3x 2 x3
L3  x   1  3 x 

2
6
3
4
2
x
x
L4  x   1  4 x  3 x 2 

3
24
5 x3 5 x 4 x5
2
L5  x   1  5 x  5 x 


3
24 120
3. 楕円関数







楕円関数と楕円積分
母数・引数の変換
ヤコビのノーム
算術幾何平均
不完全楕円積分
完全楕円積分
楕円関数・楕円積分の偏微分
楕円関数と楕円積分

難解な特殊関数


楕円関数(各種あるが)


よいライブラリが少ない
最も実用的:ヤコビの楕円関数
楕円積分
定積分=完全(complete)楕円積分
 不定積分=不完全(incomplete)楕円積分


高速ライブラリ:講師より入手可能
楕円(elliptic)関数

sn, cn, dn: ヤコビ(Jacobi)の楕円関数
sn:第1種不完全楕円積分(有理関数形)の
逆関数
x
ds
1
u
 sn x
 他は恒等式
0
2
2 2
1

s
1

k
s 


による定義

cn  u; k    1  sn  u; k  , dn  u; k   1  k sn  u; k 
2
2
引数(argument) u、母数(modulus) k
 普通は母数を省略して表記

2
楕円関数の母数


母数(modulus): k
m=k2のほうが扱いやすい


不完全楕円積分の定義式から明らか
補(complimentary)母数
k   kc  1  k

2
2

m  mc  1  m  kc
母数を陽に表記するときは表現法に注意
sn u | m  sn u; k  , etc.
楕円関数のグラフ
Jacobi Elliptic Functions: k^2=0.5
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
dn
sn
cn
0 1 2 3 4 5 6 7 8 9 1011121314151617181920
u
楕円関数の記法

グーデルマン(Guderman)記法
正弦振幅sn、余弦振幅cn、デルタ振幅dn
 正接振幅tn=sn/cn, dn(u;k)=()


グレイシャー(Glaisher)記法(分数記法)
p,q,rは(s,c,d,n)のどれか pq  u; k   pr  u; k 
qr  u; k 
 例

sn  u; k 
cn  u; k 
1
nd  u; k  
,sd  u; k  
,cs  u; k  
dn  u; k 
dn  u; k 
sn  u; k 
楕円関数の諸公式


恒等式 sn u  cn u  1, dn u  k sn u  1
特別な母数
snu,cnu,dnu  sin u,cos u,1
 k=0
2


k=1
2
2
2
snu,cnu,dnu  tanh u,sechu,sechu
特別な引数
 1

k
,
, k 
 snu, cnu, dnu   
 1 k 1 k

u=K(k)/2
 K(k)は第1種完全楕円積分

2
振幅(amplitude)関数


;k
表記  =am  u; k  逆表現 u=F  第1種不完全楕円積分(三角関数形)の
逆関数

d
1
u
0

1  k sin 
2
2
 am 
楕円関数の三角関数的表現
snu  sin  , cnu  cos  , dnu  1  k sin 
2
2
楕円関数の逆関数

逆楕円関数の基本:2変数逆楕円関数
x  rcn u; k  , y  rsn u; k   u  atn2  y, x; k 

第1種不完全楕円積分(三角関数形)による表現
atn2  y, x; k   F  atan2  y, x  ; k 

振幅関数の一般表現(下記の[]はガウス記号)
 u 
am  u; k   2 
 atan2  sn  u; k  , cn  u; k  

 4 K  k   Gauss
母数の逆数変換

母数の範囲の制限
逆数変換+虚数変換→0<k<1
 ガウス変換+ランデン変換→k~0 or k~1



k>1をk<1に
k k
sn  u; k   k sn  ku; k
1
逆数変換
1
cn  u; k   dn  ku; k
1
dn  u; k   cn  ku; k
1
1



母数の虚数変換



cn  u; k  
dn  u; k 
cn u; k
虚数変換

純虚数のkを
実数のkに
k
k k 
2
1 k
2
u  u  1 k u
2
sn  u; k  
dn  u; k  
 
dn  u; k 
1  k sn u; k
2
1
 
dn u; k
(上昇)ランデン変換

いわゆるランデン(Landen)変換


0<k<1をより大
1+ 1  k 2 sn  u ; k  cn  u ; k 
きく(=なるべく cn  u; k  
dn  u ; k 
1に近づける)
1  1+ 1  k 2 sn 2  u ; k 
2 k
sn  u; k  
k k 
dn  u ; k 

1 k
u u 
u
1 1 k
2
dn  u; k  




1  1  1  k 2 sn 2  u ; k 
dn  u ; k 
下降ランデン変換
ˆ  dn  uˆ; kˆ 
ˆ
cn
u
;
k

別名:ガウス変換
cn  u; k  
0<k<1を、より小さく
1+kˆsn  uˆ; kˆ 
k =1/2から出発して
ˆ
ˆ
ˆ
1

k
sn
u
;
k
4回の変換で k<10




sn  u; k  
ˆsn  uˆ; kˆ 


k
1+
k
k  kˆ  

 1 1 k 
1  kˆsn  uˆ; kˆ 
u
dn  u; k  
u  uˆ 
ˆsn  uˆ; kˆ 
1+
k
1  kˆ

2


2
-10
2
2
2
2
2
ヤコビのノーム(nome)


楕円関数の数値計算に便利
定義
  K  
q  exp 

 K 

完全楕円積分Kと、その補数K’
K   k   K  k 

k  1 k 2
kが大きくてもqは小さい
1

k
qe
2
0.0432
ノームの計算

|k2|<½のとき

4倍数べき展開:非常に速い収束(問:確かめよ)
q   1  2  15  150 
4


展開因子

8
k2

2 1  k   1  k 
12

2

k  1 k 2
½<k2<1のときは双対変換で補ノームへ
母数の双対変換

双対変換(補数変換 or プライム変換)
2
m  m  1  m

k  k  1 k



0<k2<1/2 → 1/2<k2<1
完全楕円積分 K , E  K , E 
ノーム
 2 
q  q  exp 

 log q 
楕円関数のノーム表現

ノームを用いた分数表現(数値計算に有利)
基本母数範囲(0<k2<1/2)に還元後を想定
 sys
k  cyc
k  yd
sn  u; k  
,cn  u; k  
,dn  u; k  
yn
yn
yn

s  sin , c  cos , 
u
2K  k 

2 q
4
k

1  2  q  q 4  q9 
1  q 2  q 6  q12 
y(q): 正規化された4つの楕円テータ関数
 問:とqをkの関数として図示せよ


平方指数展開

正規化テータ関数:収束が非常に高速
ys  1  U 3q 2  U 5 q 6  U 7 q12 
N(N+1)次:ys,yc
2
6
12
2
y

1

S
q

S
q

S
q

c
3
5
7
 N 次:yd,yn
4
9
y

1

2
T
q

T
q

T
q
2 4 6 
d
 問:高速性を
yn  1  2 T2 q  T4 q 4  T6 q 9 
体感せよ




展開係数:チェビシェフ多項式と、その変形

漸化式の活用 T  cos  n  ,U  sin n , S  Tn
n
n
n
sin
cos
楕円関数のフーリエ展開

フーリエ展開:ノーム展開ほど収束は速くない
2  q n 1/ 2
sn  u; k  
sin  2n  1 v

2 n 1
kK  k  n 0 1  q
2  q n 1/ 2
cn  u; k  
cos  2n  1 v

2 n 1
kK  k  n 0 1  q

2  q n
dn  u; k  

cos 2nv

2n
2 K  k  K  k  n 1 1  q

2q n
am  u; k   v  
sin 2nv
2n
n 1 n 1  q 

u
2K  k 
振幅関数の数値計算

振幅関数のランデン変換
vn  knun
n  sin  kn sin n 
un  vn
un1 
, vn1  un vn , n 1 
2
2
1
算術幾何平均に他ならない
 アルゴリズムの概要(問:プログラムを書け)

u0=u,v0=kuから出発し、順変換を繰り返す
 kj=vj/ujは記録しておく
 un 
1 
n1  2 tan  tanh  
 2 
 収束後はun=vn , kn=1だから

 nから出発し、kjを用いて逆変換で0まで戻す

算術幾何平均

arithmetic-geometric mean: AGM(a,b)

2変数の算術平均と幾何平均の共通収束値
an +bn
an 1 =
, bn 1  anbn
2
 一例

2

AGM 1, 1  k
 = 2K  k 
非常に収束が速い(問:確かめよ)
 平方根の多用:計算コストはかかるが…

不完全楕円積分


ヤコビ (1850):有理4次式の不定積分の基本形
ds
第1種 F  x; k   x

2
2 2
1

s
1

k
  s 
0

第2種
E  x; k   
1  s2
ds
2 2
1 k s
x
0

第3種:パラメータnを-2と表現することもある
  x; n, k   
ds
x
0
2
1

ns


2
2 2
1

s
1

k
  s 
不完全楕円積分(2)
Incomplete Elliptic Integrals: n=0.5, k^2=0.5
4
3.5
3
F(x,k)
PI(x,n,k)
2.5
2
E(x,k)
1.5
1
0.5
0
0 1 2 3 4 5 6 7 8 9 1011121314151617181920
x
不完全楕円積分(3)

三角関数を使った積分表現(ルジャンドル)


広く普及(他の表現と第1引数の違いに注意)
第1種
F  ; k   
d

1  k 2 sin 2 
0
cos2 d



第2種 E ; k   0
第3種
1  k 2 sin 2 
  ; n, k   

0
d
2
2
2
1

n
sin

1

k
sin



不完全楕円積分(4)
Incomplete Elliptic Integrals: n=0.5, k^2=0.5
24
22
20
18
16
14
12
10
8
6
4
2
0
F(phi,k)
PI(phi,n,k)
E(phi,k)
0
2
4
6
8
10
phi
12
14
16
18
20
不完全楕円積分(5)

楕円関数を使った積分表現(ヤコビ)
(再び)他の表現と第1引数の違いに注意
 他と区別するために新記号en,pnを導入




第1種 F u; k   u
第2種 E  u; k   u cn 2  v; k  dv  en u; k 
0
第3種
dv
  u; n, k   
 pn  u; n, k 
2
0 1  n sn
 v; k 
u
不完全楕円積分(6)
Jacobi Elliptic Functions: n=0.5, k^2=0.5
18
am(u,k)
16
14
12
10
pn(u,n,k)
8
en(u,k)
6
4
2
0
0
2
4
6
8
10
u
12
14
16
18
20
不完全楕円積分(7)

第2種不完全楕円積分の別表現
 K k  
E  u; k   
 u  Z  u; k 
 E k  
ヤコビのゼータ関数=E(u;k)の周期成分
Z  u; k   zn  u; k 

u
 2q sin 2v  4q 4 sin 4v  9q 9 sin 6v   
2K  k 


K  k   1  2q cos 2v  2q 4 cos 4v  2q 9 cos 6v 

不完全楕円積分(8)

特殊な場合

k=0
E ;0  F ;0  

k=1
  ; n, 0  
E ;1  sin 
tan 1

1  n tan 

1 n
 1  sin  
F  ;1  ln 

cos



1 
1
1  sin  
1
  ; n,1 
 n tan  n sin    log

n 1 
2
1  sin  

n=-k2
1
  ; k , k  
1 k 2
2

k 2 sin  cos  
 E  ; k  

2
2
1  k sin  

引数による偏微分

引数uによる偏微分
母数kとパラメータnを固定
am  u; k 
sn  u; k 
 dnu
 cnudnu
u
u
en  u; k 
cn  u; k 
2
 cn u
 snudnu
u
u
pn  u; n, k 
dn  u; k 
1
2

 k snucnu
2

u
1

n
sn
u
 u; k 

母数による偏微分
m=k2による楕円関数の偏微分(問:示せ)
sn  u | m 
am  u | m 
 cn  u | m 

m
cn  u | m 
m
dn  u | m 
 sn  u | m 
m
am  u | m 
m
sn  u | m  
am  u | m  

sn  u | m   2mcn  u | m 

m
2dn  u | m  
m

am  u | m 
 pn  u; m | m   u 
 dn  u | m  

m
2m


母数などによる偏微分

mおよびnによる偏微分(問:示せ)
en  u | m  en  u | m   u

m
2m
pn  u; n | m  1 
pn  u; m | m   pn  u; n | m  
 u 

m
2
mn

pn  u; n | m 
 pn  u; n | m   u 1  m
1


pn  u; m | m 

n
2 1  n  
n
mn



2
1  nsn  u | m  dn  u | m  
sn  u | m  cn  u | m 
不完全楕円積分の計算

一般不完全楕円積分の数値計算
(Fukushima and Ishizaki 1994b, CMDA)

G  ; nc , mc , a, b   
0

a cos 2   b sin 2 
 cos   n sin  
2
2
c
cos   mc sin 
2
2
基本公式(問:確かめよ)

 F  ; k    E  ; k   G  ;1,1  k 2 ,    ,    1  k 2 
d

 F  ; k     ; n, k   G  ;1  n,1  k 2 ,    ,  1  n    
一般楕円関数

一般不完全楕円積分の楕円関数版

基本公式
gn u; n | m, a, b  G  am u | m ;1  n,1  m, a, b 
u  en u | m  gn u;1,1 m,   ,    1 m
u  pn u; n | m  gn u;1 n,1 m,   ,  1 n   

応用その1(問:確かめよ)
K
K

zn  u | m   gn  u;1,1  m,1  ,1  m  
E
E

一般楕円関数の応用

応用その2(問:確かめよ)

注:各種の偏微分計算で必要となる
1  m  u  en  u | m   gn
m
 u;1,1  m,1, 0 
u  pn  u; n | m 
 gn  u;1  n,1  m, 0,1
n
pn  u; m | m   u
fn  u | m  
 gn  u;1  m,1  m, 0,1
m
完全楕円積分



完全積分=不完全積分の定積分
dx
第1種 K k  F 1; k  1
 
第2種

 0
E  k   E 1; k   
1
0

第3種
  n, k    1; n, k   
1  x 1  k x 
2
2
1  x2
dx
2 2
1 k x
ds
1
0
2
2
1

ns


2
2 2
1

s
1

k
  s 
完全楕円積分(2)

第3種完全楕円積分の別表現(n>0の場合)
k 2K k    
  n, k   2
 
k n  2 

n 0   ; k 
1  n   k
2
 n
  tan
1
n
k
ホイマン(Heuman)のラムダ関数
0   ; k  

E  k   K  k  F   ; k    K  k  E   ; k   

完全楕円積分(3)




楕円関数の周期:4K
 K(0)=/2, K(1)=∞
補数関係
K   k   K  k  , E  k   E  k 
k  1 k 2
母数kを省くことが多い→K,K',E,E'
ルジャンドルの関係式

 
0  , k   1  EK   EK  KK  
2
2 
完全楕円積分のグラフ
Complete Elliptic Integrals: n=0.5
4
3.5
3
2.5
2
2PI(n,k)/pi
1.5
2K(k)/pi
1
2E(k)/pi
0.5
0
0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
k
1
完全楕円積分の計算



大げさ:不完全楕円積分の計算ルーチンを呼ぶ
定番:算術幾何平均(Numerical Recipe他)
実用的:母数変換とノーム展開



例:k2<1/2の場合
2
4
第1種
 1  2q  
K  k   2 

 1 k 
第2種
 2 2   q  4q 4  9q9 
E  k   K  k   
 
4
9

K
k
1

2
q

q

q








 
母数が小さいとき

q
(母数の二乗である)mの1次までの展開
m
m
k 1
16
2
K

m
1   E
2
2
 m
2 
snu  1  cos u  sin v
4


m 2
dnu  1  sin u
2
znu

2

 m
1   u
2

 m 2 
cnu  1  sin u  cos v
 4

m
sin 2u
4
 m
enu  1   u  znu
2

4. 補間と外挿







補間とは?外挿とは?
和分、差分、差分商
ラグランジュ補間とエルミート補間
直交多項式補間
スプライン補間
平面上の補間
リチャードソンの外挿
補間とは?外挿とは?

=選点(collocation)近似:関数近似の一つ
標本点で関数値などに一致
f  xn   fn  n  0,1,...
 xn , fn   n  0,1,...

応用の対象

標本点区間の内側:補間 Interpolation
 標本点区間の外側:外挿 Extrapolation


数値微分・数値積分の諸公式の基礎
補間・外挿の分類







応用区間:内側(補間) vs 外側(外挿)
基底:多項式、三角多項式、指数多項式
基底の台:大局的 vs 局所的(スプライン)
非線形性:線形 vs 非線形(三角級数など)
微分値の利用:ラグランジュ vs エルミート
構成法:一斉的 vs 逐次的(ニュートン)
標本点の配置:等間隔 vs 非等間隔
和分

和分(mean)=平均


高階和分


f n 1/ 2  f n 1/ 2
中心和分=2者の平均  f n 
2
例:2階中心和分  2 f  f n1  2 f n  f n 1
n
用途
4
差分との組み合わせ:差分の対称化
 偏微分方程式の数値解法
 データの平滑化

差分

差分(difference)=有限な差
桁落ちに注意→差分演算子
f n  f n1  f n
 前進(forward)差分
 後退(backward)差分 f n  f n  f n 1


等間隔標本点の場合=階差
中心(central)差分
 倍幅中心差分
 2階中心差分

 fn  fn1/ 2  fn1/ 2
   fn   fn1  fn1  / 2
 2 fn  fn1  2 fn  fn1
階差表


x0 f 0
高階階差の表
前提:等間隔標本点 x1 f1 f1


x2 f 2 f 2  f 2
x3 f 3 f 3  f 3  f 3
2
3
桁落ちに注意


常微分方程式の
数値積分(多段法)
2
丸め誤差は階数の累乗で増幅される
更新アルゴリズム
f jk  k f j
f n,0 :=f n ;do(k=1,n) f n,k =f n,k-1 -f n-1,k-1
高次元の差分と和分


高次元の場合:成分別に適用
偏微分方程式の差分化に必須
中心差分、中心和分を多用
 例:2次元ラプラシアン(直交座標)の近似


仮定:等方等間隔刻み幅 x  y  h
2
2
2
2



 x y  un , m
 u  u
 2 2 
y  n ,m
 x
h2
  un 1,m  un 1,m  un ,m 1  un ,m 1  4un ,m  / h 2
差分演算子



微小だが有限な差を表現する演算子Δ
x  x  x0
定義
演算法則(問:確かめよ)
定数 c  0
 加減算   x  y    x    y 
 乗算
  xy    x  y0  x  y   x0  y    x  y

基本差分演算

問:導け
 y   y  x0  y0  x   y  x  y  x 
  

xx0
xx0
x
除算
  x 2    x  x0  x 
 2乗
3
2
2

x

x

xx

x
  
 3乗
0
0   x 
 平方根
x



 x 
立方根  3 x 
 
x  x0
x
 x    x  x    x 
3
2
3
3
0
3
0
2
基本差分演算(2)

三角関数
 x  x0   x 
  sin x   2cos 
 sin  
 2   2 
 x  x0   x 
  cos x   2sin 
 sin  
 2   2 
  tan x   1  tan x tan x0  tan  x 

逆三角関数
  tan 1 x   tan 1  x / 1  xx0  
 1  y  
1   y  x  y  x  
1   y  x0  y0  x  
  tan     tan 
  tan 

x
xx

yy
xx

yy
 

0
0
0
0




基本差分演算(3)

双曲線関数
 x  x0 
 x 
  sinh x   2cosh 
 sinh  
 2 
 2 
 x  x0 
 x 
  cosh x   2sinh 
 sinh  
 2 
 2 
  tanh x   1 tanh x tanh x0  tanh  x 

逆双曲線関数
  tanh 1 x   tanh 1  x / 1  xx0  

1  y  
1   y  x  y  x  
1   y  x0  y0  x  
  tanh     tanh 
  tanh 

x
xx

yy
xx

yy
 

0
0
0
0




基本差分演算(4)

指数関数
 x  x0 
 x 
 exp  x    2exp 
 sinh  
 2 
 2 

対数関数
 x 
  log x   log 1 

x0 


微分
 dy  d  y 
  
dx
 dx 

積分

  ydx     y  dx
基本差分演算(5)

ベクトル
成分
 加減算
 内積


 x  j    x j 
  x  y    x   y 
  x y    x x  x0  x
ベクトル積   x  y    x  y  x0   y 
基本差分演算(6)

行列
 A  jk    Ajk 
成分
T
T

A


A




 転値
 加減算   A  B   A   B
 乗算
  AB   A B  A0  B
 逆行列
1

A
1
    A   A  A
0
1
差分演算の応用

1. 万有引力の加速度(エンケ法で活用)
 1 
 1    x  x0  x  
  x 
  3      3  x   3 3  
 x0 
 r 
 r 
 r r0   r  r0  

2. ルジャンドル陪多項式の差分
P00  0, P10   sin  , P11   cos  ,
Pmm11   2m  1   cos   Pmm  cos 0  Pmm 
m
n

m

P


2n  1
n 1
m
m
m


Pn 1 

sin

P

sin


P



n
0
n 

n  m 1
n  m 1
差分商

ニュートンの差分商(Divided Difference)

定義
f n
f n  f n 1
f  xn   f  xn  , f  xn 1 , xn  

,
xn xn  xn 1
f  xn  2 , xn 1 , xn  
f  xn  k ,

, xn  
f  xn  2 , xn 1   f  xn 1 , xn 
f  xn  k ,
xn  2  xn
, ...
, xn 1   f  xn  k 1 ,
xn  k  xn
, xn 
, ...
重要:標本点列の順序に依存しない(問:示せ)
差分商(2)

差分商の性質
微分の「ある種の」拡張
 全ての標本点が同一なら微分値に比例


 f  n
f  x, , x  
 n1  n !

標本点が等間隔で順番に並んでいる場合
xn k  xn  kh  f  xn k ,
k fn
, xn  
k !h k
差分商の計算

差分商表
f ij
k
 f  xi , x j ,
, xk 
x0 F00  f 0
x1 F11  f1 F10  f10
x2 F22  f 2 F21  f 21 F20  f 210
x3 F33  f3 F32  f32 F31  f321 F30  f3210

更新アルゴリズム(桁落ちに注意)
Fn,n :=f n ;do(k=n-1,0,-1) Fn,k =  Fn,k+1 -Fn-1,k  /  x n -x k 
差分商の計算(2)

y0  x0 F00  f 0
同一点を含む場合
y1  x0 F10  f 0 F10
f k
f k
f kk  , f kkk  ,
2
6

y2  x1 F20  f1 F21 F20
更新アルゴリズム
y3  x1 F30  f1 F32 F31 F30
逆引き配列j(n)を用意
 上記の場合
yn :=x j n  ; Fn,n :=f j n  ; do(k=n-1,0,-1){
j  0   0, j 1  0
 n-k 
if  j  n  =j  k   then Fn,k :=f j n  /  n-k !
j  2   1, j  3  1




else Fn,k :=  Fn,k+1 -Fn-1,k  /  yn -y k 

ラグランジュ補間

補間の代表




多項式、関数値のみ利用、一斉的構成

ラグランジュ補間公式 PL  x  
f n Ln  x 
n
ラグランジュ多項式
x  xk
 選点直交性 Ln  xk    nk
Ln  x   
ラグランジュ補間の欠点
k  n xn  xk
等間隔補間の場合:ルンゲの現象
 一斉的構成:点数の追加に対応が困難

ルンゲ(Runge)の現象

等間隔補間の破綻の典型
標本点を増やせば増やすほど精度が劣化する
 補間される関数が実軸上に近い極を持つとき



詳細は森(2002) 、二宮(2002)を参照
問:f(x)=1/(1+x2) に、区間[0,1]において、
両端を含む等間隔標本点でのラグランジュ
補間を試み、ルンゲの現象を確認せよ

注意:(ローレンツ分布関数)f(x)の極は±i
ニュートン補間

実用的な補間


多項式、関数値のみ利用、逐次的構成
ニュートン補間公式:差分商の活用
PN  x   f  x0 
 補間点の追加
 f  x0 , x1   x  x0 
が容易
 問:ラグランジュ  f  x0 , x1 , x2   x  x0  x  x1 
補間公式との等  f  x0 , x1 , x2 , x3   x  x0  x  x1  x  x2 
価性を示せ

 f  x0 ,
n 1
, xn    x  xk  
k 0
ニュートン補間(2)

等間隔標本点の場合 xn  x0  nh
x  xn
s
h
ニュートン階差公式

s  s  1 2
 s  m  1 m
f  x  
 fn 
  f n  f n  sf n 
m 
2
m0 
 一般化二項係数
 s  m  1  s  m   s  m 



m !  s 
 m  m!

逐次多項式補間


用途:1点(x)での補間値だけ欲しいとき
エイトキン・ネビル(Aitken-Neville) の方法
更新公式(問:ニュートン補間公式から導け)
 右下端が最良近似列
x0 f 0  T00

Tn ,0  f n
Tn ,k
Tn ,k 1  Tn 1,k 1
 Tn ,k 1 
x  xn  k
1
x  xn
x1 f1  T10 T11
x2 f 2  T20 T21 T22
x3 f3  T30 T31 T32 T33
逐次有理関数補間

ブリシュ・シュテーア(Burlisch-Stoer)の方法
更新公式(パデ近似の応用)
 問:導け Tn ,0  f n

Tn ,1  Tn ,0 
Tn ,k
Tn ,0  Tn 1,0
 x  xn 1   Tn 1,0 
  1

 1 
Tn ,0 
 x  xn  
Tn ,k 1  Tn 1,k 1
 Tn ,k 1 
 x  xn  k   Tn ,k 1  Tn 1,k 1 
  1

 1 
 x  xn   Tn ,k 1  Tn 1,k  2 
エルミート補間

微分値も利用できる場合



不定積分、常微分方程式の数値解など
m ni 1
エルミート補間公式
k 
PH  x    fi Lik  x 
i 0 k 0
拡張ラグランジュ多項式
補助多項式
 定義漸化式

Li ,ni 1  x  
ik
 x
x  xi 


k!
i , n 1  x  , Lik  x  
i
k
 x  xi


j 0, j i  xi  x j
m
ni 1
ik  x   
j  k 1
 j
ik



nj
 x  Lij  x 
エルミート補間(2)

ニュートン・エルミート補間公式:実用的

例:2点3階=5次多項式
PH  x   f  x0   f  x0 , x1   x  x0 
 f  x0 , x0 , x1   x  x0  x  x1 
 f  x0 , x0 , x1 , x1   x  x0   x  x1 
2
 f  x0 , x0 , x0 , x1 , x1   x  x0   x  x1 
2
2
 f  x0 , x0 , x0 , x1 , x1 , x1   x  x0   x  x1 
3
2
エルミート補間(3)

重要:エルミート補間多項式の微分は、元
の関数の微分を、それなりに近似する

運動方程式の数値解から位置・速度の近似
多項式が構成可能:2点3階の例
PH  x   f  x0 , x1   f  x0 , x0 , x1   2 x  x0  x1 
 f  x0 , x0 , x1 , x1   x  x0  3x  x0  2 x1 
 f  x0 , x0 , x0 , x1 , x1   x  x0  x  x1  4 x  2 x0  2 x1 
 f  x0 , x0 , x0 , x1 , x1 , x1   x  x0   x  x1  5 x  2 x0  3x1 
2
エルミート補間の例
2点2階3次多項式
 両端での値と微分値の線形結合表現
P  x   f  x0  Q00    f  x1  Q10    f   x0  Q01    f   x1  Q11  



標準変数    x  x0  /  x1  x0  0    1
基底多項式(問:示せ)
Q    1     2  1 Q    1    
0
0
2
Q10     2 3  2 
1
0
2
Q11     2  1
エルミート補間の例(2)

2点2階3次基底多項式
1
0.8
Q00
Q10
0.6
0.4
0.2
Q01
0
Q11
-0.2
0
0.2
0.4
0.6
0.8
1
直交多項式補間
n

選点直交多項式

f  x    c j j  x 
j 0
定義:標本点で離散直交関係を満足
n
 w  x   x    
i 0
i
j
i
k
i
k
jk
 n1  xk   0
標本点iでの重み:wi
 実用的な標本点:高次直交多項式のゼロ点
 正規化係数 n
1 n
c j   w j j  xi  f  xi 


補間公式(問:示せ)
j
i 0
直交多項式補間(2)

代表的な直交多項式
名称
記号 補間区間 重み関数 正規化係数n
ルジャンドル Pn(x) [-1,1]
1
チェビシェフ Tn(x) [-1,1]
2/(2n+1)
ラゲル
Ln(x) [0,)
(1-x2)-1/2 n/2
(0=2, n≠0=1)
exp(-x) 1
エルミート
Hn(x) (-,)
exp(-x2)
()1/22nn!
チェビシェフ補間

チェビシェフ補間=離散チェビシェフ展開
c0 n
f  x; c     ckTk  x 
2 k 1

チェビシェフ多項式の選点直交性(問:示せ)
   
2 n
 n
 n
T j xi Tk xi   k  jk

n i 1

 n
xk
  2k  1  
 cos 

 2n 
チェビシェフ補間係数(問:示せ)
   
2 n
 n
 n
ck 
g x j Tk x j

 k n j 1
スプライン補間


スプライン(spline)=自在定規
区分多項式の一種
折れ線をn回積分→n+1次スプライン
 自然(natural)スプライン
 正規(cardinal) スプライン
 基底(basis)スプライン=B-スプライン
 3次(cubic)スプライン:実用的
 秋間(Akima)スプライン:日本の誇り

正規スプライン


選点直交性 Ck  x j    kj
欠点1:余分な決定条件が必要



N次ならN+1個:普通、両端の高階微分値を0
欠点2:計算量が多い
長所:同じ標本点→2回目以降が簡単

補間表現が直接的
f  x; x1 ,
N
, xN    f k Ck  x 
k 1
B-スプライン


基底(basis)スプライン
性質 or 定義づけ
与えられたn+1個の標本点: x0, x1, …, xn
 標本点の外側では恒等的に0
 各小区間でn次多項式
 各標本点でn-1階微分値まで連続


特徴:凸型形状
3次スプライン

特徴1:曲率最小化
3次の自然スプラインは
min
曲率のL2ノルムを最小化
 (ある意味で)最も滑らかな近似



f
x
dx







2
特徴2:そこそこ高精度&計算が容易
仮定:標本点xiで関数値と1階微分値が既知
 →[xi, xi+1]でエルミート3次補間多項式は確定
 三重対角の連立1次方程式の解法に帰着

3次スプライン(2)

3次スプラインの決定法

fkは既知なので未知数f’kを求める


標本点で2階微分が連続→ f’kの方程式


注:2階微分値を未知数として求める手法もある
問:(連立1次)方程式を具体的に導け
係数行列は三重対角→簡単な解法(後述)

注意:両端での微分値(f’0, f’N)が必要
3次スプライン(3)

重要:両端の微分値の指定法 f 0  f1  f 0
x1  x0
1.最も身近な差分商で近似
 2.近似多項式の両端での2階微分値=0



3.近似多項式の両端での3階微分値=0


3次自然スプライン
別名:Pスプライン
上記2.or 3.のいずれも条件式が追加

問:追加の条件式を導け
3次スプライン(4)

等間隔標本点の場合 (間隔=h)

連立方程式の主要部分
3
f k1  4 f k  f k1   f k 1  f k 1 
h

境界部分:自然スプライン
2 f 0  f1 

3
 f1  f 0 
h
f N 1  2 f N 
3
 f N  f N 1 
h
境界部分:Pスプライン
2
2




f 0  f1   f1  f 0 
f N 1  f N   f N  f N 1 
h
h
秋間(Akima)スプライン

特徴:そこそこの精度&計算が高速
Akima (1970, J. ACM): 3次スプラインの1種
 近傍5点のデータだけから1階微分値を推定
 →両端点では2次式に退化
 微分値推定公式:両側差分商の重み付き平均

fn   w f [ xn , xn1 ]  w f [ xn , xn1]  w  w 
w  f [ xn , xn1 ]  f [ xn1, xn2 ]
問:解析関数の補間

以下の解析関数について、補間区間をいろ
いろ変えて、下記の補間を試み、その優劣
(精度、計算時間)を論じよ




等間隔ラグランジュ、等間隔エルミート、ルジャ
ンドル、チェビシェフ、3次自然スプライン、秋間
スプライン
1
f1  x  
1. ローレンツ分布関数
2
f
x

exp

x
2. 正規分布関数
 
2 
3. 低次シュトゥンプ関数; c0,c1
1  x2
平面上の補間

安直:xy各次元ごとに補間
例:双3次スプライン、双Akimaスプライン
 難点:座標系依存性、任意領域に適用困難


本格的:デローネ(Delone)三角分割
本来:測地網の最適化理論
 双対=ボロノイ(Voronoi)多角形分割
 最低次:各三角形ごとに折れ面(1次式)補間
 次の段階:Akima 2次元スプライン bivar.f

外挿



外挿=標本点の外側で近似すること
補間は安全⇔外挿は危険
問:以下の例で確かめよ

1. sin(x)を区間[0,1]で多項式補間

標本点:両端点を含む等間隔のn+1点
2. [0,1]の外側で補間式の誤差を図示
 3. 「nを増加→誤差が急激に増加」を確認せよ

リチャードソンの外挿

非常に重要:数値計算全般に応用可能


Richardson (1910,1927)天気数値予報の始祖
有限刻み幅hの離散計算結果をh=0へ外挿
h=0:実現不可能だが、手法によらず誤差0
 hを変えて計算→計算値から補間表を構成


外挿の逐次近似アルゴリズム
多項式:エイトキン・ネビル
 有理式:ブリシュ・シュテーア

リチャードソンの外挿(2)

発想:誤差項の逐次消去
0. 仮定:計算値はhのテイラー級数で展開可能
0
p
(多項式外挿の場合)
T  h  T    T  h p 
 1. 異なるhでの計算 T  T h  T  0   T  p  h p 
 1
1
1
結果を2つ用意 T  T  h   T  0  T  p  h p 
2
2
2
 2. 主要誤差項を消去
p
p
 3. 反復
T
h

T
h
T2  T1
 0
2 1
1 2

T

h h
p
1
p
2

 T2 
 h1 / h2 
p
1

リチャードソンの外挿(3)


刻み幅変更法の王道=整数分割:h=H/n
代表的な分割数列

n
1,
2,
4,8,...,
2
,...
ロンバーグ(Romberg)列
ブリシュ列 1, 2,3, 4,6,..., 2n ,3  2n1,...
 調和(harmonic)列
1, 2,3, 4,..., n,...


hの偶数乗展開→収束が高速
関数の数値積分:ロンバーグ法
 常微分方程式の初期値問題:グラッグ法

リチャードソンの外挿(4)


T  h  T
 0

 T h2k
k 
hの偶数乗展開
k 1
ロンバーグ列、多項式補間アルゴリズム
更新公式
Tn,k 1  Tn 1,k 1
Tn,k  Tn ,k 1 
k
 数値積分で有用
4 1


調和列、多項式補間アルゴリズム
Tn ,k 1  Tn 1,k 1
更新公式
Tn ,k  Tn ,k 1 
2
 常微分方程式の
 n 

 1
数値解法で有用
nk 

5. 関数の近似

発想:標本点での一致→区間内で最も近い



距離の三大定義
最小2乗近似


「近い」=関数空間での相互距離が小さい
フーリエ展開、離散フーリエ展開
ミニマックス近似

チェビシェフ展開
関数の近似(2)

Lnノルム

距離の三大定義(n=1, 2, )
Ln  f , g   n

f  x   g  x  dx
n
L1: 通常は使われない
 L2: 最小2乗近似、計算しやすい
 L: ミニマックス近似、別名「最良近似」
→実現は難しい、チェビシェフ展開で代用

関数の近似(3)


近似基底関数の三大構成法
大局的(例:フーリエ近似)


完全局所的(例:スプライン近似)


全区間で定義、全区間で有限値
ある区間だけで定義、外部ではゼロ
部分局所的(例:ウェーブレット近似)

全区間で定義、ある区間以外でほぼゼロ
最小2乗近似

最小2乗近似の定義

近似誤差の重み付き2乗積分を最小化
min    x; p  w  x  dx   x; p   g  x   f  x; p 
p
2


w(x):重み関数=積分密度=分布関数
近似関数の決定方法
1. 「パラメータpについての偏微分=0」を解く
 2. 近似関数fがpについて線形なら
線形最小2乗法を解く

最小2乗近似(2)

例:g(x)=sin(x)を[0,1]で1次関数近似


重み=1, 近似関数 f  x   a  bx
正規方程式
A1   a   B0 
    
A2   b   B1 
 A0

 A1

解
1
1
Ak   x dx, Bk   x sin xdx
k
0
k
0
a  4  2 cos1  6sin1 ~ 0.0318
 問:導け b  6  6 cos1  12sin1 ~ 0.8558
最小2乗近似(3)
0.03
0.02
0.01
0
-0.01
-0.02
-0.03
-0.04
-0.05
0
sin(x)-(a'+b'x)
0.2
0.4
0.6
0.8
1
直交多項式近似



ランク落ちしない最小2乗近似の秘策
理想形:正規方程式の係数が対角行列
n
近似関数
f  x    c j j  x 
j 0


直交多項式
  x   x  w  x  dx   
j
k
j
近似係数公式 c  1   x  f  x  w  x  dx
j
j

j
jk
フーリエ展開


周期関数g(周期2)を最小2乗近似
近似区間:[0, 2)
2
min    ; p  d   ; p   g    f  ; p 
p

2
0
近似関数:フーリエ級数
c0 
f  ; c, s      ck cos k  sk sin k 
2 k 1
フーリエ展開(2)

三角関数の直交性(問:示せ)
2
 cos j cos k d   
k
jk
0
2
 sin j cos k d  0

2
 sin j sin k d  
jk

0
 0  2,  k  1  k  0
0

フーリエ級数の最小2乗近似係数(問:示せ)
ck 
1

2
 g   cos k d
0
sk 
1

2
 g   sin k d
0
離散フーリエ展開

周期関数g(周期2)を離散最小2乗近似

近似関数:2n+1個の三角多項式
c0 n
f  ; c, s      ck cos k  sk sin k 
2 k 1
[0, 2)中の2n個の等間隔標本点
 積分→台形則による離散和

2n


min    ; p 

p
j 0 
 n
j
2
n
j
j

n
A0  AN N 1
 Aj 
  Aj
j 0
2
j 1
N
離散フーリエ展開(2)

三角多項式の選点直交性(問:示せ)
2n
2n
i 0
i 0
 cos ji n cos ki n   k jk n  sin ji n  sin ki n    jk n
2n
 sin ji cos ki  0
 n
 n
i 0

離散フーリエ展開係数(問:示せ)
 
1 2n
 n
 n
ck   g  j cos k j
n j 0
 
1 2n
n
n
sk   g  j  sin k j 
n j 0
ミニマックス近似

ミニマックス(mini-max)近似の定義

近似誤差の最大(max)絶対値を最小(mini)化
min  max   x; p  
p

  x; p   g  x   f  x; p 
ミニマックス近似の誤差関数の特徴
1. 両端で絶対値最大
 2. 全ての極値は符号が交代で絶対値が同じ
 3. 近似関数がn次多項式→極値はn+2個

ミニマックス近似(2)


例:g(x)=sin(x)を[0,1]で1次関数近似
誤差関数のミニマックス条件 f  x   a  bx
A: 1次式なので極値は(両端を含めて)3個
 B: 両端で同じ(負)値  1    0  a  0
 C: 中間点で最大値
   x   0,   x   a


解(問:導け)
1     

a  sin   1    1  sin1 ~ 0.0300, b  sin1 ~ 0.8415
2  2   2 

ミニマックス近似(3)
0.03
0.02
0.01
sin(x)-(a+bx)
0
-0.01
-0.02
-0.03
0
0.2
0.4
0.6
0.8
1
チェビシェフ展開



フーリエ展開で変数変換 x  cos 
近似区間:[-1,1]
近似関数
c0 


チェビシェフ級数
f  x; c  
2
  ckTk  x 
k 1
チェビシェフ展開係数(問:示せ)
ck 
2
1


1
g  x  Tk  x 
1  x2
dx
6. 連立1次方程式

方程式の分類



線形1次元:自明
線形多次元:連立1次方程式の解法



線形と非線形、1次元と多次元
定番:ガウス、コレスキー、SOR、PCG
非線形1次元:複雑、ニュートン法など
非線形多次元:難解、準ニュートン法など
連立1次方程式(2)

膨大な知識の蓄積
Ax  b
A: 正方行列
3
 重要:逆行列(N の計算量)の利用を避ける


定番:良いライブラリあり
一般密行列:ガウスの消去法
 対称密行列:コレスキー法
 近似解が既知:加速緩和法(SOR)、共役傾斜法
 疎行列:形態に応じて種々の方法

右上三角行列の場合


 A11

もっとも単純
 0

後退代入

 最下行から順次に  0
A12
A22
0
bn 1  An 1,n xn
bn
xn 
, xn 1 
,
An ,n
An 1,n 1
xn  k 
1
An  k ,n  k
A1n   x1   b1 
   
A2 n   x2   b2 

   
    
Ann   xn   bn 
,
k 1


 bn  k   An  k ,n  j xn  j  ,
j 0


ガウスの消去法


時間はかかるが汎用で強力
背景思想:行列のLU分解


Lは左下三角、Uは右上三角
アルゴリズム
1. 前進消去
 2. 後退代入


Ax  b
A  LU
1
Ux  L b
1 1
xU L b
欠点:Uは直接求められない
ガウスの消去法(2)

第1段

1行目を使って
(j,1)成分を消去
 A j1 
Ajk  Ajk  
 A1k
 A11 
 A j1 
bj  b j  
 b1
 A11 
 A11

 A21


 An1
 A11

 0


 0
A12
A22
An 2
A12

A22
A11
A1n  x1   b1 
   
A2 n  x2   b2 

   

  

Ann  xn   bn 
A1n   x1   b1 
    

A2 n   x2   b2 

   
    
   xn   bn 
Ann
ガウスの消去法(3)


第2段以降:右下の小行列で繰り返す
問題点:ピボット(=対角成分)が0のとき
対策:行and/or列の入れ替え
 完全選択:絶対値最大となるように行も列も



部分選択:絶対値最大となるように行だけ


安全だがコードが複雑化
コスト・パフォーマンスに優れる←お勧め
破局:どの成分も、ほぼ0→ランク落ち

一つの処方:すべて強制的に0
LU分解(decomposition)


ガウスの消去法で同一のAを繰り返し使う
→ LとUを保存しておくと便利
LU分解の不定性→三大流儀
ドリトル(Doolittle)法 Lkk  1
 クラウト(Crout)法 U kk  1
 LDV分解 A  LDV
Lkk  Vkk  1, Djk  0  j  k 


詳細は数学ライブラリを参照

忠告:自己流はやめて既存ライブラリを使え
LU分解の応用

行列式の計算
A   LkkU kk   Dkk
k


行ピボット交換をm回した場合は(-1)m倍
逆行列の計算
ドリトル(Doolittle)法
 クラウト(Crout)法
 LDV分解法

k
QR分解

QR分解
A  QR
Q:直交行列、R:右上三角行列
T
 連立方程式の解 QRx  b  c  Q b  Rx  c
 欠点:LU分解より2~3倍遅い
 アルゴリズム:数学ライブラリを参照


応用:行列の全固有値を求める方法
QR反復
An
 次第に対角化

 Qn Rn  An1  RnQn
三重対角行列


三重対角性 Ajk  0  j  k  1
2
u
u
c 2
例:1次元拡散方程式の離散化

クランク・ニコルソン法
u j ,k 1  u j ,k
t

t
x
c  u j 1,k 1  2u j ,k 1  u j 1, k 1 u j 1, k  2u j , k  u j 1, k
 

2
2
2 

x

x
 
 




連立1次方程式の係数は三重対角  u k   u j ,k
j
→「三項方程式」という

問:A,Bを具体的に示せ
Ak 1uk 1  bk 1  Bk uk
三重対角行列(2)

ガウスの消去法


 A11

 A21
三重対角行列に応用 A   0

 0


ピボット選択を省略
 0
アルゴリズム

A12
0
0
A22
A23
0
A32
0
A33
A43
A34
A44
0
0
0
0 

0 
0 

0 


Ann 
do  k=2,n {ck-1 =1.0/A k-1,k-1;d=A k,k-1*ck-1;
A k,k -=d*A k-1,k ;bk -=d*bk-1}


x n =bn /A n,n ;do  k=n-1,1,-1 x k =ck *  b k -x k+1*A k,k+1 
コレスキー(Cholesky)法

対称行列のとき


A AAS S
T
LU分解でLとUが互いに転置行列
アルゴリズム
1. コレスキー分解
 2. 後退代入:第1段
 3. 後退代入:第2段


T
S y b
Sx  y
T
同じAを繰り返し使うとき→Sを保存
コレスキー分解


分解条件
アルゴリズム
n
S S = A   Sij Sik  Ajk
T
i 1
k 1
S11 = A11 , Skk  Akk   S 2jk
j 1
k 1
A1k
1 

S1k 
, S jk 
Ajk   Sij Sik 

S11
Skk 
i 1


破局:平方根の中が負→ランク落ち

一つの処方:強制的にSkk=Sjk=0
修正コレスキー法

LDLT分解
Aが対称→LDV分解でV=LT (問:示せ)
 アルゴリズム

j 1
k 1

1 
2
D11 = A11 , Lkj 
 Akj   Lki Dii L ji  , Dkk  Akk    Lki  Dii
D jj 
i 1
i 1


修正コレスキー法(平方根が不要)

コレスキー法でコレスキー分解→ LDLT分解
Ly  b
1
zD y
L xz
T
反復法

近似解が得やすい場合
エンケの方法:既知解が近似解
 係数行列が、ほぼ対角 or ほぼ右上三角


主な反復法:すべて1次収束
ヤコビ(Jacobi)法:非対角成分を右辺に移項
 ガウス・ザイデル(Gauss-Seidel)法



ヤコビ法の各行で、常に最新補正値を使用
加速緩和(SOR)法

ガウス・ザイデル法で修正量を最大2倍まで増加
対角優位行列


対角優位性
Akk 
n

j 1, j  k
Ajk
例:2次元ポアッソン方程式の離散化
 2u  2u
2
u

u

u

u

4
u

h
 j ,k



x
,
y
  j 1,k j 1,k j,k 1 j,k 1
j ,k
2
2
x y
連立方程式の形に書くと係数行列は対角優位
 問:A,Uなどを具体的に示せ
2

AU  h Σ
一般反復法


 n1
 n
 Mx  Nb
反復法の一般表現 x
係数行列のDEF分解 A  D  E  F
D  0  j  k
D: 対角行列
 E: 狭義の左下三角行列 E jk  0  j  k 
 F: 狭義の右上三角行列 Fjk  0  j  k 

 D11

 0
D 0


 0

jk
0
0
D22
0
0
D33
0
0
 0
0 


0 
 E21
0  E   E31




E
Dnn 
 n1
0
0
0
E32
0
0
En 2
En 3
0
 0 F12


0
0 0
0 F  0 0




0 0
0 

F13
F23
0
0
F1n 

F2 n 
F3n 


0 
ヤコビ法


一般表現 N  D , M  N  E  F 
反復公式
n
1
 k 1
xj


1

Ajj

k  
 b j   Aji xi 
i 1,i  j


Aが対角優位→ヤコビ法は収束
注意:ほとんどの場合、ヤコビ法より
ガウス・ザイデル法のほうが良い
ガウス・ザイデル法


一般表現 N   D  E 1 , M   NF
反復公式
 k 1
xj
1

Ajj
j 1
n

 k 1
k  
 b j   Aji xi   Aji xi 
i 1
i  j 1


Aが対角優位 or AとDが正値対称
→ ガウス・ザイデル法は収束
 収束が遅いとき、エイトキンの加速が有効

加速緩和(SOR)法


Successive Over Relaxation
一般表現(:加速パラメータ)
M   I  D E  1    I  D1F 
1
1
N    D  E 
1
=1 → ガウス・ザイデル法
 の決め方:様々な手法


反復公式(A,Dが正値対称、0<<2→収束)
 k 1
xj
 1
k 
 xj  
 Ajj
j 1
n

 k 1
k  
k  
  Aji xi   x j 
 b j   Aji xi
i 1
i  j 1



傾斜(gradient)法

最小2乗問題へ置換 min    r


別名:逐次最小化法
 k 1
解の修正 x


k 
多様な決定法(後述)
直線上の最小探索

k 
 x  k p
探索方向決定

x
解表現(問:示せ)
r  b  Ax
   , 
p   F r 
k
2
j
, p
j
k 
j
k 
p r
min    k   k 
k 

p Ap
傾斜法(2)

探索方向の代表的な決定方法

最急降下(steepest descent)法 p  r


k 
k 
Aの固有値が違うと収束が非常に遅い
共役(conjugate)傾斜(CG)法
探索方向ベクトル:グラム・シュミットの直交化
 定番:前処理付共役傾斜(preconditioned CG)法



準(quasi)ニュートン法
詳細は「非線形最適化問題」で
共役傾斜法



共役(conjugate)=直交(orthogonal)
k 
 j
直交条件 p Ap  0  j  k 
グラム・シュミットの直交化(問:示せ)
k 
p


r
k 
 k 1
 k 1p
 k 1
k 
p
Ar
 k 1   k 1
 k 1
p
Ap
有限(N次元ならN)回探索で解に到達
欠点:丸め誤差に弱い⇒一度は没落
前処理付共役傾斜法

PCG (=Preconditioned CG)法


背景:クリロフ部分空間法(森 2002)
発想:Aを前処理+共役傾斜法
Ax  b  By  c c  C b B  C A  C
1

1

T 1
y  CT x
前処理の思想 B  I
B=Iなら1回の探索で解に到達(問:確かめよ)
 例:Aが対称→不完全コレスキー分解


不完全コレスキー分解共役傾斜 (ICCG)法
不完全コレスキー分解


=不完全LDLT分解
LDLT分解の欠点


発想:Aij=0なら強制的にLij=0


Aが疎でもLが密(問:示せ)
アルゴリズムの詳細→数学ライブラリを参照
Aが非対称→不完全LDV分解

詳細は福井(1999)を参照
7. 非線形方程式


見かけの複雑さ
非線形方程式の一般論





グラフによる近似解法
低次代数方程式
1次元非線形方程式
高次代数方程式
連立非線形方程式
見かけの複雑さ


重要:見かけに惑わされない
実は単純な例
3sin x  cos 2 x  1  0
変数変換で簡易化 s  sin x
 問:解け(落とし穴に注意)


真に難しい例
x  1  sin x
非線形方程式の一般論

解の存在と唯一性
多重解の困難
 近似解の重要性:特にグラフ解法の有効性
 解の囲い込み:既存知識、グラフ解法、多分法



反復解法→収束判定の難しさ
解の妥当性の吟味
物理的に正しい解か?
 得られた解は求めていた解か?

多重解の例

例:x-1=9cos(x)
10
5
0
-5
-10
-15
-10
-5
0
5
10
グラフによる近似解法

適当なグラフ・ソフトを用意


お勧め:GNUPLOT
手法
1. 解の探索:グラフを描き、近似解を読み取る
 2. 収束判定
 3. 精密化:描画範囲を狭める


例題 x  1  sin x
グラフ解法の例
10
x-1
5
0
sin(x)
-5
-10
-15
-10
-5
0
5
10
グラフ解法の例(2)
1.2
sin(x)
1
0.8
x-1
0.6
0.4
0.2
0
1
1.2
1.4
1.6
1.8
2
グラフ解法の例(3)
1
0.98
x-1
0.96
0.94
sin(x)
0.92
0.9
0.88
1.9
1.92
1.94
1.96
1.98
2
方程式解の収束判定

条件1: 方程式の値(本筋) fn   f


絶対判定基準(f0は代表値)  f
 f0 
条件2: 解の変化量(要注意) xn1  xn   x

三大判定基準
絶対:x0は代表値  x  x0 
 相対:解が大変動する場合  x  xn 
 複合  x  max  x0 , xn  


判定微小量の設定例   10 100

1次元非線形方程式








準1次方程式
2次方程式
高次代数方程式
2分法と多分法
逆補間
セカント法とミューラー法
ニュートン法とハレー法
デュラン・ケルナー法
準1次方程式

ほとんど1次方程式 ax  b    x 
例:ケプラー方程式でeが小さいとき
E  M  e sin E
 単純反復法が有効
1
 1次収束
x0  a b
 連立方程式でも適用可
1


より良い方法
xn 1  a b    xn  
微小項の計算が面倒→セカント法
 微小項が微分可能→ニュートン法

2次方程式


落とし穴:根の公式
2
b  b  4ac
x

2
困難1:近接根(4ac~b )
2a


やむをえない精度劣化
c
x2 
ax1
方法1:根と係数の関係を活用
困難2:準1次方程式(ac~0)




ax  bx  c  0
2
第1根←分子の絶対値が大きくなる符号の方
方法2:分子の有理化
実用的:ニュートン法
x
2c
b
b  4ac
2
近接2根の困難

一般の非線形方程式でも同様




近接2根の近傍では、ほぼ2次方程式だから
x 
 根の相対誤差:計算精度の半分の桁数
現象1: 根の精度劣化 x
現象2: 丸め誤差→複素根となる危険性
現象3: 反復法の収束次数の低下

例:ニュートン法で2次→1次
2次方程式(2)

問:なるべく正確な実根を求めよ
真に2次の場合 x 2  3x  1  0
 準1次方程式の場合

0.000001x  2 x  1  0
2

近接根の場合
x  2.000001x  1  0
2

解の吟味が必要な例(変数変換せよ)
4sin x  2 cos 2 x  1
高次代数方程式

根の公式:限定的かつ非実用的


3次:カルダノ、4次:フェラリ (問:調査せよ)
実用的な方法
単根:ニュートン法、ハレー法
 実係数方程式の2複素共役根:ベアストウ法
 全根

安直:(単根探索+単項式の除算)などの繰り返し
 本格的:一斉に求める→連立法(DK法など、後述)

多項式の四則演算



式の演算=係数ベクトルの演算
加減算:自明
乗算:単純
N
M
N M



i
j
k
  ai x     b j x    ck x
 i 0
  j 0
 k 0

係数公式
k
ck   ai bk i
i 0
多項式の四則演算(2)

除算:割り切れない限り不可能
割り切れるとき(その1):多項式÷単項式
N 1
 N

i
j
a
x

x

x

b
x

0
j
 i  
j 0
 i 0

 係数の漸化式:乗算公式の逆(問:確かめよ)
b N-1:=a N

do(i=N-2,0,-1) bi :=a i+1 +bi+1*x 0 
多項式の四則演算(3)

割り切れるとき(その2):多項式÷2次式

2次式の根が実根でないときに有効

i
2
j
A
x

x

rx

s

B
x
  j
 i  
j 0
 i 0

 係数の漸化式:乗算公式の逆(問:確かめよ)
BN-2 :=A N ;BN-3:=A N-1 +r*BN-2
N
N 2
do  i=N-4,0,-1Bi :=Ai+2 +r*Bi+1 +s*Bi+2 
2分(bisection)法

遅いが確実 xL  x  xR f L : f  xL  , f R : f  xR 


仮定:囲い込みが終了
アルゴリズム
1)中点で評価
 2)収束判定
 3)囲い込み更新


1次収束
xn1  x ~ xn  x / 2
xM   xL  xR  / 2
x M :=0.5*  x L +x R  ;f M :=f  x M 
収束判定
if  f L *f M <0 
then x R :=x M ;f R :=f M 
else x L :=x M ;f L :=f M 
多分(multisection)法

グラフ解と等価:2分法より遅いが、より確実



並列計算機ならメリットあり
主な用途:複数根の分離
アルゴリズム:多点評価+囲い込み更新
 枝分かれ do  k=1,n-1 x k :=  k*  x 0 +x n   /n;f k :=f  x k 
に注意
do  k=0,n-1{
if  f k *f k+1 <0 

then m++;x 0m :=x k ;f1m :=f k+1
逆補間(regula falsi)



2分法の変形:確実に収束 x  xL f R  xR f L
N
fR  fL
反復公式:1次補間式を解く
x N :=  x L *f R -x R *f L  /  f R -f L 
 問:求めよ
f N :=f  x N 
アルゴリズム


2分法とほぼ同様
1次収束(問:示せ)
 f  
xN  x 
xL  x  xR  x 


2f 
収束判定
if  f L *f N <0 
then x R :=x N ;f R :=f N 
else x L :=x N ;f L :=f N 
エイトキン(Aitken)加速


用途:1次収束を超1次収束に加速
2
xn 1  xn 
 注:3点が必要

xn 1  xn 1 
xn 1  2 xn  xn 1
加速公式
 元の数列が1次収束 xn1  x   C   n  xn  x 

収束条件を連立させて解く


xn  xn 1  C  xn 1  xn 1 
微小項nを無視
加速後は超1次収束 xn1  x

問:示せ
xn 1  xn 1  C  xn  xn 1 
xn1  x

 n   n1
C  C  1
2
セカント(secant)法

逆補間の変形:高速だが不安定




ニュートン法より高速だが収束の保証が無い
反復公式:逆補間と同じ
アルゴリズム:囲い込みの確認を放棄
x n+1:=  x n *f n-1 -x n-1*f n  /  f n-1 -f n 
超1次収束
f n+1:=f  x n+1 
 次数~1.618
 f  
収束判定
xn1  x 
x

x
x

x



n
n 1

2f 
n++
セカント法の加速

エイトキン加速の応用


注:4点が必要
加速公式
xn 1  xn 1 
 xn 1  xn 
2
xn 1  2 xn  xn  2
元の数列の収束 xn1  x  C   n  xn  x  xn1  x 
 連立させて解く
xn 1  xn 1  C  xn  xn 1  xn 1  xn 1 
 微小項nを無視
xn  xn 1  C  xn 1  xn 1  xn  2  xn 1 
 問:加速後の収束条件式を求めよ

ミューラー(Muller)法

セカント法の拡張:複素根に強い
反復公式:2次ニュートン補間式を解く
 2次方程式の解:絶対値が小さい方を選ぶ


複合は分母の絶対値
xn1  xn 
が大となる方
bn
2cn
b  4ancn
2
n
an  f  xn2 , xn1, xn , bn  f  xn1, xn   an  xn  xn1  , cn  fn


アルゴリズム:セカント法と同様(問:導け)
超1次収束:次数~1.84
ニュートン(Newton)法





高速だが不安定:安定な初期値が必要
微分値が容易に計算可能なら有効
fn
反復公式
xn1  xn 
f n
 1次テイラー展開を解く
アルゴリズム
2次収束(問:示せ)
 f n 
2
xn1  x 
  xn  x 
 2 f n 
x n+1:=x n -f n /f n
 :=f   x n+1 
f n+1:=f  x n+1  ;f n+1
収束判定
n++
ニュートン法(2)


有効(=微分値が容易に計算可能)な例
 ケプラー方程式 f  E   E  e sin E  M  0
安定な初期値
区間内で下に凸→上界、上に凸→下界
 例:ケプラー方程式の場合
 前処理:modと鏡像変換 0  E    E  
0


2次収束の例外

多重根の場合:1次収束(問:示せ)
ニュートン法(3)

等価な最適化問題 f  x   0  opt f N  x 
ニュートン写像関数
f  x
 仮定:2階微分が0とならない f N  x   x 
f  x
 問:等価であることを示せ


安定な初期値=部分集合での最適値

例:ケプラー方程式に対するNijenhuisの初期値
M e 
 M
E0  min f N  E   min 
, M  e,

  
1 e
1 e 

E 0, , 
2


ニュートン法(4)

応用1:平方根 f  x   x  a  0
2
反復公式(問:示せ)
1
a 
xn 1   xn  
2
xn 
 1<a<4の場合に
安定な初期値の例 x  min  1  a , 3  a ,1  a 
0






部分集合(1,3/2,2)
2
2 4
3
4
問:aが上記以外の場合は、どうすればよいか
問:n乗根のニュートン反復公式と安定な初
期値の例を示せ。
ニュートン法(5)

応用2:平方根の逆数
f  y   y2  a  0
yn 1  0.5 yn  3  ayn2 
反復公式(問:示せ)
 1<a<4の場合に安定な初期値の例

部分集合
  3   a   9   27a   3   a  
y0  min       ,    
,     
  4   16   8   128   2   2  
(1/2,3/4,1)

平方根の逆数yが求まれば平方根自体は
a*yで求まる
 除算が不要なので、より実用的

ニュートン法(6)
応用3:逆数 f  x   x  a  0
1


除算が低速(orハードで計算不能)なら有効

好例:逆行列の計算
xn1  xn  2  axn 
反復公式(問:示せ)
Xn1  2Xn  Xn AXn
 反復公式:行列の場合
 問:行列Aが対角優位のとき、対角成分Dの逆
行列から出発して、ニュートン法で逆行列A-1を
求めるプログラムを書き、他の方法と比較せよ

ニュートン法の加速

エイトキン加速の応用その2


3点が必要
xn1  xn1 
2  xn1  xn 
 xn1  xn1 
2
3
 3  xn1  xn 
加速公式の導出
xn 1  x   C   n  xn  x 
元の数列の収束
 連立させる
xn1  xn1  C  xn  xn1 
2
微小項nを無視
xn  xn1  C  xn1  xn1 
2



2
得られる2次方程式をニュートン法で解く
2
ニュートン法の変形

多重根の場合の収束速度の低下防止策

変形方程式を導入
すべて単根(問:示せ)
 桁落ちに注意


f  x
0  g  x 
f  x
ニュートン公式を適用:どんな場合も2次収束

ハレー法(後述)
との類似性に注意
 g0
 f0
x 

f 0f 0
g0
f 0 
f 0
ベアストウ(Bairstow)法

発想:「代数方程式を2次式で割った剰余=
0」の方程式をニュートン法で解く
q  x  b x
p  x   q  x   x 2  rx  s   Ax  B  0
n2
k 0


未知数:2次式 x2-rx-s の2係数(r,s)
1
方程式
 A A 


A=B=0
ニュートン法
r
 r   r
     
 s n 1  s n  B

 r
n  2k
k

s  A 
  
B   B n

s n
ベアストウ法(2)

A,Bを求める漸化式(問:求めよ)
b0  a0 ; b1  b0 r  a1
bk  bk 2 q  bk 1r  ak (k  2,
, n  2)
A  bn3q  bn2 r  an1; B  bn2 r  ai

偏微分を求める漸化式(問:上記から導け)
b1,r  b0 ; b2,r  b1,r r  b2
b2,q  b0 ; b3,q  b2,q r  b1
, n  2) bk ,q  bk  2,q q  bk 1,q r  bk  2 (k  4, , n  2)
A
B
A
B

b
q

b
r

b
;
 bn  2,q q  bn  2
 bn3,r q  bn2,r r  bn2 ;
 bn2,r q
n 3, q
n  2, q
n 3
q
q
r
r
bk ,r  bk 2,r q  bk 1,r r  bk 1 (k  3,
ハレー(Halley)法





E. Halley:周期彗星の予言で有名
高速で、かなり安定(詳細は不明だが)
2階までの微分値が計算可能なら最強
fn
巧妙な反復公式
xn 1  xn 
f n f n
f n 
3次収束(問:示せ)
xn1  x
 2 f n f n 3  f n2 
3

  xn  x 
2
 12  f n 

2 f n
ハレー法(2)

公式の導出法:2次テイラー展開
f 0
2
0  f 0  f 0x   x  
2

形式解の右辺にニュートン公式を代入
 f0
 f0
x 

f 0
f 0f 0
f 0  x f 0 
2
2 f 0
ハレー法(3)


応用1:平方根 x 2  a
 xn2  3a 
xn 1   2
 xn
 3 xn  a 
応用2:立方根 x3  a
 xn3  2a 
xn 1   3
 xn
 2 xn  a 

応用3:放物線軌道のケプラー方程式
3
t
t M
3
3M  6Mt  t  2t  t
tn1 
3  6Mtn  3tn2  2tn3
2
n
3
n
4
n
5
n
デュラン・ケルナー
(Durand-Kerner)法

高次代数方程式


N
k 0
j 1
発想:f(zj)=0を連立して解く


全根を一斉に求める
f  z    ak z k  aN   z  z j 
N
ニュートン法→2次法、ハレー法→3次法
ニュートン反復公式のDK近似
 n 1
zk
 n
 zk




 

z 
 
f  z 
a   z   z  
f zk
n
n
n
k
f zk
k
n
0
k
n
n
j
j k

問:同様にして3次法の公式を導出せよ
問1: ケプラー方程式

問:様々な(M,e)に対して、各種の方法で楕
円ケプラー方程式を解き、14桁正しい解を
得る計算時間を比較せよ
f  E   E  e sin E  M  0

参考
0  M  , 0  e 1
M e 
 M
M  E  min 
, M  e,

1

e
1

e


解Eの存在範囲
 単純反復法は、大きなe(>0.6)では遅い収束
 sinとcosの同時計算法を使え

問2: ルジャンドル
多項式の根


問:様々な方法で、n次ルジャンドル多項式
の根を求めよ(13次までの解は既出)
参考
Pn+1の根はPnの隣りあう根で挟まれる
 漸近形
2

1

Pn  cos   
sin   n    
n sin 
2
4

→正根の近似解  n
  4k  1  

xk  cos  
 
  4n  2  
問3: ルジャンドル
多項式の極値点


問:様々な方法で、Pnの極値点= Qnの根を
求めよ(14次までの解は既出)
参考
Qnの根はPnの隣りあう根で挟まれる
 漸近形
1
2n

1

Qn  cos  
cos   n    
sin  
2
4

→正根の近似解  n
  4k  3  

xk  cos  
 
  4n  2  
連立非線形方程式




一般論:解の存在、唯一性、適合性
主な解法
ニュートン法:2次収束
準(quasi)ニュートン法:超1次収束
=多次元セカント法
 定番:逆行列版ブロイデン公式


次元の呪い:囲い込みの困難
ニュートン法(7)



多次元:連立1次方程式の解が必要

ニュートン方程式 J  x  x  f  x 
n
n
n

ヤコビ行列(Jacobian)
f
J x  x
x
計算量は多いが2次収束
暴走(overshooting)の危険性
ニュートン法(8)

ニュートン法の暴走防止策


減速(damping)係数 xn1  xn  n xn
係数決定法
直線探索 min f  xn  xn   n
 大局的最適解の発見は困難
 一つの方法   m   
 m  1, 2,
2


初期値発見の困難

連続変形法




仮定1:方程式の分解 f  x  g  x   h  x   0
仮定2:基本方程式の解が既知 g  y   0  y
方程式を少しずつ変形して解を逐次求める
f  ; z   g  z   h  z   0   0  1; z  y  x
計算量は多いが確実

問:ケプラー方程式(e=0.9)を連続変形法で解け
f  e; E    E  M   e sin E  0 e  0  0.9; E  M  ?
準ニュートン法

ヤコビ行列の計算を避ける Kn xn  fn
K: ヤコビ行列の近似行列
fn  f  xn1   f  xn 
 セカント条件 K n 1xn  f n
 注意:セカント条件だけでは一意に決まらない


ブロイデン(Broyden)の更新公式
K n 1  K n

f n  K n x n   x n


x n
2
初期値:何も情報がなければ単位行列 K 0  1
準ニュートン法(2)

ブロイデン更新公式の難点



連立方程式を解く必要性
発想:ヤコビ行列の逆を近似 xn
逆行列版ブロイデン更新公式
 Lnfn
 x n  L n f n   x n  L n
L n 1  L n 
x n  L n f n 

初期値:何も情報がなければ単位行列 L0  1
8. 最適化問題


目的関数Fを最大化/最小化 opt F  x
最適化問題の2大分類
制約条件なし:一般、簡易
 制約条件あり:難解


注意:最大化問題→最小化問題
max F  x   min   F  x  

以降では最小化問題として扱う
1次元最適化問題






1次元非線形方程式の解法と深い関連
グラフによる解の探索
2分法と多分法
黄金分割法
放物線近似(ブレント法)
セカント法と接線交差法
非線形方程式との関係

前提:関数の最適値→(偏)微分のゼロ点




注1: 関数が(一部区間でも)微分不可能なら不成立
注2: 微分可能でも逆は不成立(問:反例を挙げよ)
注3: 最大か最小かの判定は常に必要
囲い込み(bracketing)の基本
 方程式の解:逆符号対  xL , xR : xL  xR , F  xL  F  xR   0

最適化問題:関数値の凸三つ組(convex triplet)
 xL , xM , xR  : xL  xM
 xR ,  F  xL   F  xM    F  xR   F  xM    0
1次元最適化問題(2)

最適解の存在と唯一性
多峰性の困難
 近似解の重要性:特にグラフ解法の有効性
 解の囲い込み:既存知識、グラフ解法、多分法



桁落ちの困難:微分値が計算不能な場合
最適解の妥当性の吟味

物理的に正しい解か?求めていた解か?
グラフによる解の探索

探索の手順
1. まず十分広い範囲でグラフを描く
 2. グラフから最適解の近似存在範囲を読
み取り、描画範囲を狭める
 3. 単峰性が確認されるまで2.を繰り返す
 4. 単峰性を確認したら数値探索へ




例題: min 0.01 x 2  x  cos3x 


グラフ探索の例
2
1.5
0.01(x^2-x)+cos(3x)
1
0.5
0
-0.5
-1
-10
-5
0
5
10
グラフ探索の例(2)
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0
0.01(x^2-x)+cos(3x)
0.5
1
1.5
2
グラフ探索の例(3)
-0.986
-0.988
-0.99
0.01(x^2-x)+cos(3x)
-0.992
-0.994
-0.996
-0.998
-1
1
1.02
1.04
1.06
1.08
1.1
2分法

遅いが確実


xL  x  xR
仮定:囲い込みが終了 FL : F  xL  , FR : F  xR 
アルゴリズム
x M :=0.5*  x L +x R  ; FM :=F  x M 
1)中点で評価 収束判定
 2)収束判定
if  FL <FM or FL <FR 
 3)凸三つ組更新
then x R :=x M ; FR :=FM 


1次収束
xn1  x ~ xn  x / 2
else x L :=x M ; FL :=FM 
多分法


グラフ探索と等価:2分法より遅い、より確実
アルゴリズム:最小値判定+凸三つ組更新

1、2回実施した後は2分法か数値探索へ
FL =F  x L  ; k M :=0; x M :=x L ; FM :=FL ; h:=  x R -x L  /n
do  k=1,n  x k :=x L  k*h; Fk :=F  x k  ;
if  Fk <FM k M :=k; x M :=x k ; FM :=Fk 
if  k M =0  then x R :=x L +h; FR :=F1
elseif  k M =n  then x L :=x R -h; FL :=Fn-1
elseif  Fk M 1 <Fk M +1  then x L :=x M -h; x R :=x M ; FL :=Fk M 1 ; FR :=FM 
else x R :=x M +h; x L :=x M ; FR :=Fk M 1 ; FL :=FM 
最適化の収束判定





凸性の確認:最重要  Fn1  Fn  Fn1  Fn2   0
条件1: 関数の変化量 F  F  
条件2: 解の変化量 x  x  
判定基準:絶対値、相対値、複合判定
判定の閾(しきい)値  F  Fn 
n1
n1

桁落ちへの配慮
n
n
F
x
2 Fn
2 F
x 

Fn
Fn

最適化の収束判定(2)

微分値が解析的に計算可能な場合

もちろん凸性は確認
 Fn1  Fn  Fn1  Fn2   0

条件3: 微分値の値 F   
n1
F

判定の閾値の例  F   F0 

添字0は代表値の意味
黄金分割法

2分法の変形
xn1  1  r  xn  rxn1
Fn  Fn1
仮定:単調減少列
 両端を比較して、より解に近い方に近づく
 近接法:黄金分割(golden section)比
r  5  1 / 2  0.61803 x M := 1.0-r  *x n  r*x n-1; FM :=F  x M 




収束判定
アルゴリズム
if  F <F 
1次収束(問:示せ) then x :=x
M

xn1  x
1 r  xn  x
n
n+1
M
; Fn+1:=FM 
else x n+1:=x n ; Fn+1:=Fn ; x n :=x M ; Fn :=FM 
ブレント(Brent)法

凸三つ組のデータを基に、放物線近似

近似式の頂点を逆に解く(問:求めよ)
xn  xn 2   Fn  Fn 1    xn  xn 1   Fn  Fn 2 

xn 1  xn 
2  xn  xn 2  Fn  Fn 1    xn  xn 1  Fn  Fn 2  
2

2
アルゴリズムの核心
(xn+1, Fn+1)の値による凸三つ組の更新
 問:自ら導け(参考:Numerical Recipe 10.2)

セカント法と接線交差法



微分値が容易に計算できる場合
凸三つ組→凸境界(=微分値の逆符号対)
 xL , xR  ; FLFR  0
セカント法
微分値=0をセカント法で解く
 問:凸境界の更新アルゴリズムを示せ


接線交差法
xM
FR xR  FLxL    FR  FL 


FR  FL
関数値も利用
 更新点:凸境界の両端からの接線の交点

多次元最適化問題


間接法:非線型方程式
F
opt F  x  
0
の求解問題に変換
x
直接法(最適値を直接探索)の構成
初期値・初期探索方向の決定
 直線探索=1次元最適化問題に帰着
 収束判定
 探索方向dの更新:多種多様な方式

直接法の分類

探索方向の更新方法による分類
最急降下(steepest descent)法
 ニュートン法
 共役(conjugate)方向法: Powell
 共役傾斜法: Fletcher and Reeves (1964)
 準ニュートン法



別名:セカント法、可変計量(variable metric)法
定番:逆行列版BFGS公式
最急降下法

Fを最適値のまわりに1次テイラー展開
F  x  F  xˆ   g  xˆ  x 


傾斜(gradient)ベクトル: g
F
g
x
x  x  xˆ
最急降下法 (Cauchy 1848)
各点でgの計算が必要
dn  g  xn 
 探索方向=最急降下方向
 非常に遅い(1次収束)が確実(大域的収束)

ニュートン法

2次テイラー展開
F  x  F  xˆ   g  xˆ  x  x  H  xˆ  x  / 2 


ヘッセ(Hesse)行列H:対称
ニュートン法
2 F
H jk 
x j xk
各点でgとHの計算が必要→膨大な計算量
 探索方向=ニュートン方向 H  x  d  g  x 
n
n
n
 連立1次方程式を解く必要→コレスキー法
 速い(2次収束)が非減少探索方向の可能性

準ニュートン法(3)

ニュートン法の欠点
「探索方向=減少方向」が非保証
 Hの計算量が膨大←例外:最小2乗法では緩和


発想:ヘッセ行列を正定値行列Dで近似


高速性+大域的収束性
代表的なヘッセ行列の更新公式
Davidon-Fletcher-Powell (DFP)
 Broyden-Fletcher-Goldfarb-Shanno (BFGS)

準ニュートン法(4)
Dndn  gn

準定番:BFGS更新公式

定番:逆行列版BFGS更新公式
Dn xn    Dn xn  gn  gn

Dn1  Dn 

g n xn
 Dn xn  xn
ヘッセ行列の逆行列の近似行列E dn  Engn
 連立1次方程式を解く必要なし
En g n   xn  xn   En g n   g n  En g n   xn  xn

E 
 1

En1
n
g n xn


g n xn

 g n x n
制約付き最適化問題

2種類の制約条件


不等式 Gk  x  0 等式 H j  x   0
F, Gk, Hj が全て線形→線形計画法
Linear Programming (LP)
 定番:シンプレックス(単体)法


それ以外=非線形計画法(NLP)
等式制約:ラグランジュの未定乗数法
 不等式制約:KKT条件付未定乗数法

未定乗数法

ラグランジュ関数 L  x,    F  x     j H j  x 
j


ラグランジュの未定乗数(multiplier): 
最適条件
 L

L
 optF  x  , H j  x   0    x  0,   0 
j



結論:制約なし最適化問題に帰着
未定乗数法(2)

拡張ラグランジュ関数
L  x, ,    F  x    k Gk  x     j H j  x 
k

j
Karush-Kuhn-Tucker (KKT) の5条件
L
L
L
 0,
 0,
 0, k  0, k Gk  x   0
x
k
 j

結論:乗数条件付きの最適化問題に帰着

詳細は鈴木他(1994)を参照
9. 最小2乗法

最もよく使われる最適化問題


Gauss (1801) 失踪小惑星セレスの軌道推定
目的関数

2
1
F  p    wm  g  tm ; p   g m 
2 m
差の2乗和
パラメータ:p, 近似関数:g, 座標:t, データ:g
 一般には近似関数もデータもベクトル値
1
 重み:w ← データの誤差
wm  2
m
 このとき「カイ2乗和」

線形最小2乗法

線形:近似関数gがパラメータcについて線形
座標の依存性は不問
 g=基底関数の線形和
g  t; c    ck k  t 


k
基底関数:座標1次元の例
n
多項式
 フーリエ級数
g  t ; c    ck t

k
k 0
n
g  t; s, c     sk sin k t  ck cos k t 
k 0
見かけの非線形性


注意:基底関数(パラメータ)を選ぶとき、
見かけの複雑さに惑わされないこと
いくつかの例
1. 三角関数の振幅と位相→加法定理で展開
Asin t     S sin t  C cos t
 2. 2次式の標準形→多項式形

a  x  x0   d  ax 2  bx  c
2
線形最小2乗法(2)
n

座標1次元の例(続き) g t; c    c
k 0
k
exp  k t 
n
指数関数
g  t; s, c     sk e  t sin k t  ck e   t cos k t 
 減衰振動
k 0
n
 各種の直交関数
g  t; c    ck X k  t 

k
k
チェビシェフ多項式:|t|<1
k 0
 ルジャンドル多項式:有限区間(|t|<1)
 ラゲル多項式:半無限区間(0<t)
 エルミート多項式:無限区間


上記の組合せ
線形最小2乗法(3)

m
座標2次元の例
n
g  x, y; c    cnk x k y nk
n 0 k 0
2元多項式

 n
g   ,    J n     Cnm cos m  Snm sin m 
 円盤調和関数
 a  m0
n 0


座標3次元の例
3元多項式
 球面調和関数


a
g  r, ,      
n 0  r 
m
n
k
g  x, y, z; c    cnkj x k y j k z nk
n 0 k 0 j 0
n 1 n
m
P
 n sin  Cnm cos m  Snm sin m 
m0
観測方程式

理論値と観測データの等式
N
 c  t   g
k
k
m
k 1
一般には過剰方程式系
=未知数Nより方程式の数Mが(大幅に)過剰


観測方程式
理論や観測がベクトル値→1列に並べ直す
 係数行列はN行M列(=縦長の長方形)


略記法
Ac  g
Akj   j  tk 
m
正規方程式


最小化:目的関数のパラメータ偏微分=0
 n

F M
  wm   c j j  tm   g m  k  tm   0
ck m1  j 0

正規方程式:pについての連立1次方程式
Bc  b

M
Bkj   wm j  tm  k  tm 
m 1
BはN行N列正方実対称行列

Bkj  B jk
標準解法:コレスキー法
M
bk   wn k  tm  g m
m 1
分散共分散行列 C   A A 
T


分散共分散(variance-covariance) 行列
(安直な)パラメータ推定誤差 pk   Ckk
  2Fmin / W
標準偏差
W   wm
 仮定:パラメータ間の分離が良
m
 =共分散行列の逆行列が対角行列


正しい推定誤差:Cの固有値が必要
誤差楕円:Fの2次形式表現=高次項を無視
 なかなか実行されない←コスト高

1
正規方程式(2)

注意:正規方程式のランク落ち
不適切な基底関数の選択
 不必要に多数の基底


対策1:基底関数の分離(直交性)の向上


選点直交多項式(例:チェビシェフ多項式)
対策2:無理やり解く
(特にムーア・ペンローズの)一般逆行列
 観測方程式の係数行列を特異値分解

一般(generalized)逆行列
1 
A A

 0
 定義:Aが縦長→A+は横長

観測方程式の解 Ac  g  c  Ag
+
 注意: A は一意でない


ムーア・ペンローズ(Moore-Penrose)のA+

以下の4条件を満たす一般逆行列

AA A  A,



A AA  A ,
 AA
 T
 AA  ,

A A

T
 A A
特異値分解(singular
value decomposition)

特異値分解
A  USV
U U 1
T
T
U:N行M列、S:M行M列対角、V:M行M列直交
 アルゴリズム:数値計算ライブラリを参照


ムーア・ペンローズの一般逆行列の構築

 T
A  VS U Skk  1/ Skk  Skk  0 
 Skk  0 

0
問:4条件を確かめよ
 分散共分散行列

 A A   V S
T
1

 2
VT
平均と標準偏差


最も単純な最小2乗法
スカラー・データに対する定数値の推定

答は重み付き平均値(問:確かめよ)
G
cˆ 
W

M
W   wm
m 1
M
G   wm g m
m 1
標準偏差の計算(データ数は十分大と仮定)
2Fmin

W
2
2 Fmin
G
 G2 
W
M
G2   wm g m2
m 1
直線回帰

1次関数(1次元)の当てはめ g t; c   c0  c1t
例1:測光データから背景光(スカイ)を除去
 例2:スペクトル・データから連続光部分を除去
 正規方程式
M
M
 B0 B1   c0   b0 
k
k

      Bk   wmtm , bk   wm g mtm
m 1
m 1
 B1 B2   c1   b1 
 解
b0 B2  b1B1
b1B0  b0 B1

cˆ0 
B0 B2  B
2
1
, cˆ1 
B0 B2  B12
平面回帰

1次関数(2次元)の当てはめ
g  x, y; c   c0  cx x  cy y
例:CCD画像データから背景光部分を除去
 正規方程式
 B0 Bx By   c0   b0  B   w , B   w x , B   w y ,

    B  w x ,B  w x y ,B  w y ,



 Bx Bxx Bxy   cx    bx 
B B
  c y   by  b   w g , b   w g x , b   w g y
B
xy
yy  
  
 y

M
0
m 1
M
m
x
M
M
xx
m 1
m 1
y
m
m 1
M
2
m m
xy
M
0
m m
m 1
M
m 1
m m
m
yy
M
m
m
x
m 1
m
m 1
m
2
m
m
m
M
m
m m
y
m 1
m
例外値(outlier)の除去


観測データの場合は非常に重要
除去判定条件
g  tm ; c   g m  3
3シグマ(3)管理
 閾(しきい)値(例:負値) g  t ; c   G
m
min


アルゴリズム上の工夫


データの採用・不採用を示す論理配列を用意
重要:例外値を除去 → cやが変化

新たな除去データがなくなるまで繰り返す
データの部分的集約

膨大なデータ:正確に計算→計算コスト大


例:能動的観測値(探査機などとの通信)
正規点(normal point)
=座標が近いデータ群を集約する代表点
 最も簡単な例:均質データの場合

正規点=N個のデータ群の座標の重心
 代表値=N個のデータ群の値の平均
 重み=個々のデータの重みのN倍

データの部分的集約(2)

実用的な方法
1. 座標について一定範囲ごとにデータを分離
 2. 分離後の各データ群に対し、低次多項式で
最小2乗法(例外値の除去も含む)による当て
はめを実施
 3. 正規点=除去後の重み付き重心
 4. 代表値=正規点における推定多項式の値
 5. 重み=除去後のデータ群の重みの総和

データの部分的集約(3)

例:(等間隔・等重み)5点データを2次近似
間隔:h, データ: (g-2 , g-1 , g0 , g1 , g2)
2
 2次近似式 g  n   a  bn  cn
 最小2乗法→係数と代表値(問:求めよ)

2 g 2  g 1  2 g 0  g1  2 g 2
2 g 2  g 1  g1  2 g 2
b
,c 
10
14
3g 2  12 g 1  17 g 0  12 g1  3g 2
a
35
三角回帰

三角関数(1次元)の当てはめ
g t; S , C   S sin t  C cos t
例:時系列データから特定周波数成分を検出
 正規方程式
 BSS BSC   S   bS  B   w sin t , B   w cos t ,

    
 BSC BCC   C   bC  B   w sin t cos t ,
 注意
B  B  W   w b   w g sin t , b   w g cos t

M
M
2
SS
m
m 1
2
m
CC
m 1
m
m
M
SC
M
SS
CC
m 1
m
m 1
m
m
M
m
S
m 1
M
m
m
m
C
m 1
m
m
m
ペリオドグラム
(Periodogram)


ロム(Romb) 太陽の黒点データの解析
発想:三角回帰で解(S,C)をの関数と解釈
2
2
 パワー・スペクトル P    S    C  


位相スペクトル     atan2  S   , C  
特徴
長所1:非等間隔データでもOK → 天文向き
 長所2:高い分解能
 欠点:計算量は多い

拡張ペリオドグラム

発想:基底関数を三角関数以外に拡張
Harada and Fukushima (2003)
 非線形パラメータpを固定して、線形パラメータc
について線形最小2乗法を解く

n
g  t; c, p    ck k  t , p 
k 1

パワー・スペクトル=目的関数の減少値の2倍
P  p   2F  2  F  ck  cˆk   F  ck  0  
拡張ペリオドグラム(2)

(指数関数的)単調減衰

指数関数+定数→減衰率と定常値を推定
g t; c   c0  c1 exp  t 
例:超新星の減光曲線
M
 正規方程式
B 

 B0

 B1
0
w
m
M
, Bk   wm exp  k tm ,
m 1
m 1
B1   c0   b0 
M
M

   
wm g m , b1   wm g m exp  tm 
B2   c1   b1  b0  
m 1
m 1
拡張ペリオドグラム(3)

減衰振動 g t; S, C    S sin t  C cos t  exp  t 
2次元(減衰率と周波数)スペクトルの例
B   w exp  2 t  sin t ,
 正規方程式
 BSS BSC   S   bS  B   w exp  2t  cos t ,

    
 BSC BCC   C   bC 

M
2
SS
m 1
m
m
m
m
m
M
2
CC
m 1
m
M

問:1次項A+Btを加えた
ときの正規方程式を
求めよ
BSC   wm exp  2 tm  sin tm cos tm ,
m 1
M
bS   wm g m exp   tm  sin tm ,
m 1
M
bC   wm g m exp  tm  cos tm
m 1
拡張ペリオドグラム(4)

ケプラー運動 g t; A, S, C   A  S sin E t   C cos E t 
2次元スペクトルの別の例
 例:伴星(惑星?)による視線速度変動
 離心近点角Eはケプラー方程式を解いて求める

E  e sin E  nt  E t 


非線形パラメータは2つ:平均運動nと離心率e
問:正規方程式の表現を求めよ
非線形最小2乗法



2種類の非線形最小2乗法:広義と狭義
狭義:全てのパラメータについて非線形
広義:一部のパラメータについて非線形


例:近似関数が非線形基底の線形和の場合
F  F  c, p  g t; c, p   c  t; p 
一例:フーリエ級数(周波数も未知)+1次項
n
g  t; a, b, s, c, ω   a  bt    sk sin k t  ck cos k t 
k 0
広義から狭義へ

広義最適化問題の簡易化


=広義から狭義への還元
発想:非線形パラメータpを固定
解釈:cの線形最小2乗解はpの関数 cˆ  cˆ  p 
 同一視:Fはpだけの関数

Fˆ  p   F  cˆ  p  , p 

結論:狭義の非線形最小2乗法に帰着
広義から狭義へ (2)

傾斜ベクトル

ヘッセ行列
Fˆ  F 


p  p  c  cˆ
 2 F 
  2 F   cˆi 
 2 Fˆ

  2 
 

p pm  p pm c cˆ

c

p

p
i 
i
m  c  cˆ 

 2 F
 

i , j  ci c j

問:上記表現を示せ
  cˆi   cˆ j 
 



p

p
 m 
c cˆ 
広義から狭義へ (3)


(非線形パラメータpに関する)線形最小2
乗解の係数cの偏微分が必要
正規方程式が解ける場合:簡単
←正規方程式と類似の連立1次方程式の解
M
 c   b   B 
Bkj
 j
wm
B     
 tm ; p  k tm ; p 
 c pi  2
pi
m1

p

p

p
    

M

b
 k
k
 問:特異値分解が必要な
  wn
 tm ; p  g m
p j m1
p j
場合の手法を導け

非線形最小2乗法(2)



以降は狭義非線形最小2乗法に限定
二ュートン法 F  F0  g0 p  p  H0p / 2
ニュートン探索方向d H pn  dn  g pn 
傾斜ベクトル F
 g 
  wm 

 ヘッセ行列
pk
m
 pk 

 g
2 F
H jk 
  wm 
p j pk
m
 p j
 g t
m
; p   gm 
m
  g    2 g
 
  
m  pk m  p j pk


  g  tm   g m 

m
非線形最小2乗法(3)


ニュートン法の欠点:計算量が大きい
対策1:ガウス近似(H計算の簡易化)
 g   g 
Hの表現中で
H jk   wm 
 


第2項を無視
m
 p j m  pk m
 注意:残差が小さくないと危険


対策2:ハートレー(Hartley)法
減速係数の導入 pn1  pn  ndn
 直線探索の必要性:コスト高

非線形最小2乗法(4)

対策3:マルカート(Marquardt)法

ガウス近似と単位行列の混合で近似→安定化
 g
H jk   j jk   wm 
 p
m
 j

さまざまな混合係数
レーベンベルグ
(Levenberg)係数
 ナッシュ(Nash)係数

  g 
 

m  pk m
 g
   wm 
 p
m
 j
L
j
 Nj  1   Lj
  g
 
m  p j


m
10. 数値微分







等間隔標本点の場合
非等間隔標本点の場合
補間多項式の微分
近似関数の微分
1階微分と高階微分
桁落ちの影響
数値微分の活用法
等間隔数値微分


等間隔標本点の場合(数値積分結果など)
階差表現の導出:演算子法の活用
微分演算子 Dy  f
 シフト演算子 Syn  yn1
 テイラー展開と後退差分の演算子解釈

hD 

h k 
yn 1   yn  S  
 ehD
k!
k 0 k !
k 0

k

k
yn  yn  yn1   1 S 1  1 exp  hD
等間隔数値微分(2)

微分の差分・和分表現(問:示せ)

後退差分と前進差分
hD   log 1   log 1   


中心差分と和分
1
1
hD  2sinh  / 2  2cosh 
倍幅中心差分と中心差分和分
2
1
1   
hD  sinh    
sinh  
2
1  2 / 4
等間隔数値微分(3)

微分の中心差分和分表現(問:確かめよ)
2

 n !
2
sinh     


2
n  0  2n  1 !
1  / 4
hD 

1   
2

2 n

3 5 7
     

6 30 140

1階微分の中心差分和分公式(問:示せ)
2次
 4次
 6次
 8次

2hfn  fn1  fn1
12hfn    fn2  fn2   8  fn1  fn1 
60hfn   fn3  fn3   9  fn2  fn2   45  fn1  fn1 
840hf n    f n  4  f n  4   32  f n 3  f n 3   168  f n  2  f n  2 
+672  f n 1  f n 1 



等間隔数値微分(4)

2階微分の中心差分表現(問:確かめよ)
2
4
6
8
10
12










h2 D 2  4 sinh 1      2   



12 90 560 3150 16632
 2 


2階微分の中心差分公式(問:示せ)
h2 fn   fn1  fn1   2 fn
2次
2
 4次 12h fn    fn2  fn2  16  fn1  fn1   30 f n
2
 6次 180h fn   f n3  fn3   21 fn2  fn2   255  f n1  f n1   470 f n
 f   64  f  f   784  f  f 
 8次 5040h f     f

2
n
n4
n4
n 3
 7616  f n 1  f n 1   13790 f n
n 3
n2
n2
非等間隔数値微分

非等間隔標本点の場合(点xnでの微分値)
ニュートン補間公式を微分してx=xnを代入
 実用的ヒント:標本点はxnから近い順に交互に

f n  f  xn , xn 1   f  xn , xn 1 , xn 1   xn  xn 1 
 f  xn , xn 1 , xn 1 , xn  2   xn  xn 1  xn  xn 1 
 f  xn , xn 1 , xn 1 , xn  2 , xn  2   xn  xn  2  xn  xn 1  xn  xn 1 

 f  xn , xn 1 , xn 1 ,
m 1
, xn  m    xn  xn  k  xn  xn k 
k 1
 f  xn , xn 1 , xn 1 ,
m 1
, xn  m , xn m   xn  xn  m    xn  xn  k  xn  xn k  
k 1
桁落ちの影響

数値微分の原理的困難
似たような値の差→桁落ちが深刻
→刻み幅を小さく取れない→刻み幅の最適値
 対策:近似誤差~丸め誤差←刻み幅を調節


重要:n次公式の計算精度はn/(n+1)どまり

例:1階微分の
2次中心差分
fn 
f n 1  f n 1  1 3 2

f n  

  fn h 
2h
2h

 6
hopt  3 3 f n  2 / f n  ,  f n min  3 3 f n
3
2
f n   2
3
数値微分の活用法

微分公式の検査

複雑な関数の常微分・偏微分
→公式の導出やコード化を間違いやすい


例:陰関数、超高次元、多数回の変数変換
検査目的:倍精度2次公式(誤差12桁)で十分
 f 
  
 x  y ,



f x 1  3  , y,
 

 f x 1  3  , y,
2x
 
3
注意:事前スケーリングが重要

数値微分の活用法(2)

常微分・偏微分の代用
超高精度は不要→計算コストは大きくない
 常微分



2次公式の評価~関数評価2回
偏微分
発想:原点値を共有
 1次公式の成分毎コスト=関数評価N/(N+1)回




f x 1   , y,  f  x, y,
 f 
  
x 
 x  y ,

11. 数値積分

線積分
有限区間:ニュートン・コーツ、ガウス型
 等間隔データ:中心補間型公式
 周期関数の1周期積分→台形則が最良
 無限・半無限区間:ガウス型、二重指数関数型
 特異点の処理→二重指数関数型


定番:外挿法(ロンバーグ法)
有限区間の数値積分

既知関数の定積分を線形結合で近似
b
N
a
k 1
I   f  x  w  x  dx   Ak f  xk 


まずは区分化
標本点の分布
M
I 
am
 f  x  w  x  dx
m 1 am1
等間隔:ニュートン・コーツ(Newton-Cotes)公式
 直交多項式の根:ガウス型公式

ニュートン・コーツ公式


等間隔の多項式補間公式を全域で積分
最も重要:台形則
h


2点公式
I2 
f  a   f  a 

2
0
1
他の公式(あまり使われない)
3点:シンプソン1/3
 4点:シンプソン3/8

h
I 3   f  a0   4 f  a1/ 2   f  a1  
3
h
I 4   f  a0   3 f  a1/ 3   3 f  a2 / 3   f  a1  
8
ガウス型公式

直交多項式による補間公式を積分
N
I N   wk f  xk 
k 1

評価点=高次直交多項式の根
ルジャンドル多項式:P13まで既に記述
 チェビシェフ多項式:三角関数表現


N 1
2
1
1
重み係数
   j  xk  
 問:導け w
k 0 k
k
ガウス型公式(2)

直交多項式による分類
Pn: ルジャンドル、有限区間(端点値含まず)
 Pn: ロバット、有限区間(端点値含む)
 Tn: チェビシェフ、周期関数に有効
 Ln: ラゲル、半無限区間積分
 Hn: エルミート、無限区間積分

ガウス・ルジャンドル公式

有限区間、評価点は区間の内部点

一次変換で基本区間[-1,1]に変換後
1
N
1
k 1
 
I N   f  x  dx   wk f xk



N
N
 
 
評価点xk:PNの根
N  2 

2 1  xk


n




次数:2N
wk 
2
N
  

nP
x
重みwk
 N 1 k 
ロバット(Lobatto)型公式

有限区間、評価点は両端点と内部点
一次変換で基本区間[-1,1]に変換後
1
2  f  1  f  1 N  N 
N 
  f  x  dx 
  wk f xk
 N  2  N  1 k 1
1

I N 2



 
評価点xk:PN+1 = QN+1の根
2  N  1
次数:2N+1 w N  
k
2
N


重みwk
 N  2  PN  xk 
中心補間型公式

有限区間、評価点は両端点と等間隔外部点
I 2 N 2 
xm1
N
 f  x  dx   w
k 0
xm


N 
k
 f  xm  k   f  xm 1 k  
次数:2N+3
積分の中心差分和分表現(問:導け)

hD


2sinh 1  / 2  1   2 / 4
  2 11 4 191 6 2497 8 14797 10
92427157 12
  1  





 12 720 60480 3628800 95800320 2615348736000



中心補間型公式(2)

中心差分和分公式(問:導け)
3次=台形則 I2   h / 2 fn1  fn 
I   h / 24  13  f  f    f  f  
 5次
I   h /1440  802  f  f   93  f  f   11 f  f  
 7次
I   h /120960  68323  f  f   9531 f  f 
 9次

4
6
8
n 1
n2
n
n 1
n 1
n2
n
n 1
n
n 1
n 3
n2
n2
n 1
1879  f n3  f n2   191 f n 4  f n3 

11次
I10   h / 7257600   4134338  f n1  f n   641776  f n 2  f n1 
162680  f n3  f n2   28939  f n 4  f n3   2497  f n5  f n4 
数値的畳み込み


L  f  s   F  s   
1
c  i
c  i

多項式→微分
L
分数式→積分
L

f  t  exp  st  dt
ラプラス変換
1
L
F
t

f
t

ラプラス逆変換      2 i 
0



1
1
F  s  exp  st  ds
df
 sF  s   f  0    dt
t
 F s 

  0 f   d
 s 
畳み込み L  FG t    f * g t    f t  g t   d
有理式との数値的畳み込み→数値微積分
1

0
数値的畳み込み(2)





 f    F     f t  exp  it  dt
フーリエ変換
1
F
F
t

f
t

フーリエ逆変換      2 
F
1



多項式→微分

分数式→積分
F
F
1
1
F   exp  it  d 
df
 i F     dt
 F   
t

   f   exp i  t     d
 i     

畳み込み F  FG t    f * g t    f   g t   d
応用例:Shirai and Fukushima (2000, AJ)
1

周期関数の積分
2
基本:重み1での1周期積分 I  f  x  dx
0
 オイラー・マクローリン展開
ba
B h




 f  x  dx  h  f  a  kh     2 j !  f b   f  a  h  n

n
b
k 0
a



2j
2j
2 j 1
2 j 1
j 1
周期関数→上記展開の誤差項は非常に小
結論:台形則は非常に高精度
I  h f  a  kh 
周期性より、和が簡単になる
 応用:離散フーリエ変換公式→FFT

n 1
k 0
周期関数の積分(2)

例:ベッセル関数の積分定義
J 2m  y  
1

cos  2mx  y sin x  dx


0

台形則による近似(問:示せ)
1  n 1
 2mk
 k
J 2 m  y   1   cos 
 y sin 
n  k 1
 n
 n

 
 
 
問:標本点数n=4,8,16として、適当なm, yで
積分の精度を調べよ
二重指数関数型公式

Takahashi and Mori (1974)


無限区間のオイラー・マクローリン展開→被積
分関数が無限遠で急減少なら台形則は最良
発想:無限区間に変数変換+台形則
I   f  x  dx  I 


N

n  N
 f   y  dy  h  f   nh    nh 
最適な変換=変換後が2重指数関数的
f   y      y   exp  a exp y 
二重指数関数型公式(2)


端点が特異点のときは効果絶大
有限区間の場合
準備:1次変換で積分区間を[-1,1]に
 変数変換


  y   tanh  sinh y 
 積分公式
2


I
h
2
N

n  N

cosh nh


f  tanh  sinh nh  
2
  cosh 2   sinh nh 



2

二重指数関数型公式(3)

半無限区間の場合
積分区間 [0,]
 変数変換


  y   exp  sinh y 
2

 積分公式

I

h
2
N

n  N





f  exp  sinh nh   exp  sinh nh  cosh nh
2

2


問:右の積分を求めよ I 


0
dx
1
 x
3
二重指数関数型公式(4)

半無限区間の場合(その2)

被積分関数が既に指数関数的に減少する場合
f  x   g  x  exp  x 
変数変換


  y   exp   y  exp   y 
 積分公式
2


I
h
2

N

n  N


  

f  exp  nh  exp  nh    nh  exp  nh  1  exp  nh 
2
 2


問:右の積分を求めよ

I

0
exp   x 
1 x
dx
二重指数関数型公式(5)

無限区間の場合
積分区間 [-,]
 変数変換


  y   sinh  sinh y 
2

 積分公式

I

h
2
N

n  N





f  sinh  sinh nh   cosh  sinh nh  cosh nh
2

2


問:右の積分を求めよ I 




dx
1 x
2

3
ロンバーグ(Romberg)法

発想:台形則の誤差はhの偶数累乗展開
0
2
4
6
I  h   I    I   h2  I   h4  I   h6 

台形則の積分結果をh→0へ多項式外挿
I n ,k  I n ,k 1   I n ,k 1  I n 1,k 1  /  4  1
k
分割数は倍々で増やす(ロンバーグ列)
 試験積分の計算に前回の結果を利用可能


積分点数が可変、コードが簡単→定番
ロンバーグ法(2)

第一段
I1,0

ba 

f a  2 f

4 

 a b 

  f  b 
 2 

第二段以降:前段の結果を再利用
h1  0.5b  a 
 hを半減 hn1  0.5hn

→新たな標本点の総和の追加ですむ
2n
I n1,0  I n,0  hn1Sn1 , Sn1   f  a   2 j  1 hn1 
j 1

1階差分は直接計算可能 In1,0  In,0  hn1Sn1
問:数値積分

さまざまな方法で次の定積分の値を求めよ

パラメータの値は許容範囲内で適宜与えよ

I1 
1

x



2 5/ 2
1
I 2   1  x
dx
I3  e  



cos x
1  e cos x
dx
1



2 3/ 2


3
dx
n
x
I 4  c, n   
dx
x
1  ce
0
12. 常微分方程式の
数値解法 I. 1段法

差分法の分類
1段法と多段法
 陽公式と陰公式




外挿法
特殊な2階常微分方程式専用の方法
適応的積分法=パラメータの可変化

刻み幅、次数、外挿の段数、…
常微分方程式



特別な場合:線形、変数分離形など
一般形と初期条件 dy
 f  y , t  y  t0   y0
dt
高階→連立
dn y
 f  y, y,
n
dt

 dyn
,t   
 f  y, t  ,
 dt
y  y, y,

, y
n
y ,y ,
特殊な2階

運動方程式に頻出
0
1
2
dyk
,
 yk 1 ,
dt
, yn 
d x
 f  x, t 
2
dt
dy0

,
 y1 
dt

1次元線形の場合

一般形

解の積分表現


dy
 a t  y  b t 
dt
y t0   y0
t


1
y t  
 y0   b   z   d 
合成積の登場
z  t  

t0
補助解の表現
t

z  t   exp   a   d 
 t0

1次元線形の場合(2)


定数係数:連立定数係数の場合の基本解
dy
 ay  b  t 
y t0   y0
dt
解の積分表現:合成積=畳み込み

数値的畳み込みが可能
y  t   y0e
 a  t t0 
t
  b   e
t0
a  t  
d
連立線形の場合

一般の場合

変分方程式などで頻出
 非線形の場合と同様に解くしかない
定数係数行列の場合 A t   A

dy
 A t  y  b t 
dt
1
1)係数行列の対角化
PAP  D
z  Py, c  Pb
 2)変数変換
 3)変換後、成分別に解く dz
 Dz  c  t 
dt

変数分離形

1次元、自律的


変数の主従を逆転→数値積分に帰着
dy
dt
1
dy
 f  y 

t  
dt
dy f  y 
f  y
1次元、積の形

媒介変数を導入→数値積分+補間に帰着
dy
ds
ds
1
dy
 f  y  g t    g t  , 
 s   g  t  dt  
dt
dt
dy f  y 
f  y
オイラー(Euler)法


Euler (1768)
陽(explicit)公式
yn1  yn  hf n
dy y
別名:前進オイラー

h
 思想:微分を1階差分商で近似 dt
2
 特性:1段、1回、1点、1次  y

n1  h f n
 局所誤差:fの微分値


略記法 yn  y tn  , fn  f  yn , tn  , h  tn1  tn
差分法の特性指数

段数(step): n


回数(stage): r


fを何回計算するか=主コスト
点数(point)


未来予測に必要な(現在及び過去の)総ステップ数
いくつ異なる時点(t)での情報が必要か=副コスト
次数(order): p


公式がhの何乗まで正しいか
誤差定数C:誤差の比例定数
 yn
yn
 Ch
p 1
差分法の特性量

線形安定領域(Linear Stability Region)
線形問題 y’=y に対して「安定な」
複素パラメータ h の領域
 安定:成長率(=|yn+1/yn|)が1以下


周期安定境界(Interval of Periodicity)
特に単振動問題 y’=iy において安定な
最大位相HMAX  h
 =線形安定領域の虚軸上の限界

オイラー法(2)

陽公式の線形安定条件
yn1
yn1  yn  h yn 
 1   1
yn


安定領域=中心-1/2,半径1/2の円の内側
周期安定境界


虚軸上に安定領域なし→ HMAX=0
結論:陽公式は非常に不安定
オイラー法(3)

陰(implicit)公式
yn1  yn  hf n1
別名:後退オイラー
 硬い常微分方程式で多用
 陰公式=直接には計算できない

一般:予測子修正子法(後述)で収束するまで解く
 連立線形:付随(adjoint)連立1次方程式を解く

dy
 A t  y  b t 
dt
1  hAn1  yn1  yn  hbn1
硬い常微分方程式


硬い(stiff)=大きさが異なる固有成長率
(or減衰率)が複数ある方程式系
 f  x, y  
 簡単な例 d  x   x


    8    
dt  y  10 y 
 g  x, y  
振動的に(oscillatory)硬い=大きさが異
なる固有振動数が複数ある方程式系

典型例:太陽系天体の軌道運動

イオ:1日、海王星:約300年
オイラー法(4)

陰公式の線形安定条件
yn1
1
yn1  yn  h yn1 

1
yn
1 


安定領域=中心1/2,半径1/2の円の外側
周期安定境界


虚軸上は全て安定→ HMAX=
結論:陰公式は非常に安定
テイラー級数法

オイラー法の第一の拡張
高階微分の計算が容易でない限り非効率
積分公式
h 2  df  h 3  d 2f 

高階微分の計算例


y n 1  y n  hfn 

    2 
2  dt n 6  dt n
df  f   f 
  f   
一般には非常に複雑 dt  y   t 
  2f   f   df    2f 
d 2f   2f 

 ff  2 
f       2 
2
dt
 yy 
 yt   y   dt   t 
ルンゲ・クッタ法


Runge(1895), Heun(1900), Kutta(1901)
オイラー法の拡張の第二の方向
思想:中間点(=内部点)で試行錯誤
 長所1:事前の準備が不要=自動出発可能
 長所2:安定領域が広い
 短所1:コスト高
 短所2:高次公式が複雑かつ導出が困難
 短所3:陰公式は非実用的(例外あり)

ルンゲ・クッタ法(2)

ルンゲ・クッタ法の特徴


1段、多数の内部点、多数回評価
低次の公式
2次2回2点陽公式:中点則 (Runge 1895)
 2次2回1点陽公式:修正台形則
 2次2回1点陰公式:台形則
 3次3回3点陽公式:Heun (1900)
 3次3回2点陽公式:Kutta (1901)
 4次4回2点陽公式:The ルンゲ・クッタ

低次ルンゲ・クッタ法


2次2回2点陽公式:中点(midpoint)則
h
tn p  tn  ph
yn 1/ 2  yn  f n
2
fn1/ 2  f  yn1/ 2 , tn1/ 2 
yn 1  yn  hf n 1/ 2
2次2回1点陽公式:修正台形則(ホイン法)
yn 1  yn  hf n

h
yn 1  yn 
f n  f n 1
2

f n1  f  yn1 , tn1 
低次ルンゲ・クッタ法(2)

2次2回1点陰公式:台形(trapezoidal)則

「硬い」方程式で重宝、対称公式
h
yn 1  yn   f n  f n 1 
2

fn1  f  yn1, tn1 
連立線形の場合:付随連立1次方程式を解く
h
h
 h



 1  A n1  y n1   1  A n  y n   bn  b n 1 
2 
2
 2


 3次以上の陰公式は非実用的
低次ルンゲ・クッタ法(3)

3次3回3点陽公式:ホインの3次公式
k1  hf n
k1
z2  yn  ; k2  hf  z2 , tn 1/ 3 
3
2k 2
z3  yn 
; k3  hf  z3 , tn  2 / 3 
3
1
yn 1  yn   k1  3k3 
4
低次ルンゲ・クッタ法(4)

3次3回2点陽公式:クッタの3次公式
k1  hf n
k1
z2  yn  ; k2  hf  z2 , tn 1/ 2 
2
z3  yn  k1  2k2 ; k3  hf  z3 , tn 1 
1
yn 1  yn   k1  4k3  k3 
6
低次ルンゲ・クッタ法(5)

4次4回2点陽公式:Theルンゲ・クッタ法


RK4と略記
k1  hf n
特徴:2点公式
k1
z2  yn  ; k2  hf  z2 , tn 1/ 2 
2
 f の計算で t 依存
k2
部分が重い場合は、 z3  yn  ; k3  hf  z3 , tn 1/ 2 
2
依存部分の再利用 z4  yn  k3 ; k4  hf  z4 , tn 1 
により約2倍高速化
1
yn 1  yn   k1  2k2  2k3  k4 
6
低次ルンゲ・クッタ法(6)

RK4:連立線形の場合

注目:An+1とbn+1は次のステップで再利用可能
k1  h  A n y n  b n 
1
z 2  y n  k 1 ; k 2  h  A n 1/ 2 z 2  b n 1/ 2 
2
1
z 3  y n  k 2 ; k 3  h  A n 1/ 2 z 3  b n 1/ 2 
2
z 4  y n  k 3 ; k 4  h  A n 1z 4  b n 1 
1
1
y n 1  y n   k 1  k 4    k 2  k 3 
6
3
ルンゲ・クッタ法(3)

誤差(陽公式):ほぼp+1次のテイラー展開項
理由:p次陽公式の係数は、テイラー展開のp次
までの項を再現するように決定するから
5
 例:RK4、単振動  y
 h 
n


線形安定領域

yn

120
p次(p<5)陽公式では、以下の領域
zn
 1 z 

n 0 n
p
zp

1
p
ルンゲ・クッタ法(4)

陽公式の安定領域: RK1~4 (複素上半面)
3
RK4
2.5
RK3
2
1.5
RK2
1
Euler
0.5
0
-3 -2.5 -2 -1.5 -1 -0.5 0 0.5
ルンゲ・クッタ法(5)

陽公式の周期安定境界


決定方程式(p<5)
周期安定境界の例
1  iH MAX 
iH MAX 


Euler: 完全不安定
 RK2: 完全不安定
HMAX  3 1.732
 RK3
 RK4
HMAX  2 2 2.828


高次になるほど周期安定境界は大
p
p
1
ルンゲ・クッタ法(6)

台形則(2次陰公式)の線形安定条件
yn1 1   / 2
h
yn 1  yn 

1
 yn1  yn  
2
yn
1  / 2


安定領域=複素左半面=実部が0または負
周期安定境界


虚軸上は全て弱安定→ HMAX=、摂動に弱い
結論:台形則は、かなり安定
高次ルンゲ・クッタ法

ブッチャー(Butcher)の次数障壁


5次6回、6次7回、7次9回、8次11回、…
主な高次公式
Fehlberg: RKF4(5), RKF7(8)
 Verner: DVERK


Dormand and Prince: DOPRI5, DOP853
DOP853の係数は複雑
 http://www.unige.ch/~hairer/software.html

高次ルンゲ・クッタ法(2)


一般表現
J回N(M)次
N:計算子
 M:比較子
k1  hf n
 2  tn  c2 h; z2  yn  a21k1 ; k2  hf  z2 , 2 
 3  tn  c3h; z3  yn  a31k1  a32 k2 ; k3  hf  z3 , 3 


比較子
最終段公式
 刻み幅制御
に用いる

j 1
 j  tn  c j h; z j  yn   a ji ki ; k j  hf  z j , j 
i 1
J
yn 1  yn   b j k j
j 1
J
yˆ n 1  yn   bˆ j k j
j 1
高次ルンゲ・クッタ法(3)

一般表現の係数表とRK4の例

注:RK4には比較子は無い
c2 , a21
c3 , a31 , a32
cm , am1 , am 2 ,
am,m 1
y,
bm
bˆ
yˆ ,
b1 ,
bˆ ,
1
b2 ,
bˆ ,
2
m
1 1
,
2 2
1 1
,0,
2 2
1, 0, 0,1
1 1 1 1
y, , , ,
6 3 3 6
刻み幅(step size)制御

適応(adaptive)積分法:自動積分の一種


可変刻み幅、可変次数、可変次数可変刻み幅
最も簡単な適応積分法:刻み幅制御
仮定1:局所誤差が測定可能(例:比較子)
p 1
 仮定2:局所誤差の次数が既知  y  h
 制御法:変更後の局所誤差増加率~制御値

1/ p
 ch 
 c  h  h 

h
 y 
y
刻み幅制御(2)

アルゴリズムの一例
ˆ
r:=(1.0+abs(y))*c*h/(abs(y-y)+EPS)
if (r > 2.0 || r < 0.5)
then {h:=h*exp(log(r)/p)}
注1:複合誤差を採用
 注2:頻繁な変更を避けている点に注目


安全策

変更後の誤差が制御範囲内に入ることを確認
1 1
,
4 4
3 3 9
, ,
8 32 32
12 1932 7200 7296
Fehlberg (1969)
,
,
,
13 2197 2197 2197
4(5)次6回6点
439
3680 845
1,
,  8,
,
216
513 4104
 比較子のために、
一度、中点(1/2) 1 , 8 ,2, 3544 , 1859 , 11
まで戻ってk6を評 2 27 2565 4104 40
価する点に注意 y, 25 ,0, 1408 , 2197 , 1 ,0
216 2565 4104 5
比較子のほうが
16
6656 28561 9 2
yˆ ,
,0,
,
, ,
高次→損
135 12825 56430 50 55
RKF4(5)



DOPRI5



Dormand and
Prince (1980)
5(4)次6回6点
巧妙な方法

k7は次ステップ
のk1として再利
用可能
1 1
,
5 5
3 3 9
, ,
10 40 40
4 44 56 32
, ,
,
5 45 15 9
8 19372 25360 644448 212
,
,
,
,
9 6561 2187
6561 729
9017 355 46732 49 5103
1,
,
,
,
,
3168 33 5247 176 18656
35
500 125 2187 11
1,
,0,
,
,
,
384 1113 192 6784 84
35
500 125 2187 11
y,
,0,
,
,
, ,0
384 1113 192 6784 84
5179
7571 393 92097 187 1
yˆ ,
,0,
,
,
,
,
57600 16695 640 339200 2100 40
DVERK


Verner (1978)
6(5)次8回7点
1 1
,
6 6
4 4 16
, ,
15 75 75
2 5 8 5
, , ,
3 6 3 2
5 165 55 425 85
,
, ,
,
6 64 6 64 96
12
4015 11 88
1, ,  8,
,
,
5
612 36 255
1 8263 124 643 81 2484
,
,
,
,
,
,0
15 15000 75 680 250 10625
3501 300 297275 319 24068
3850
1,
,
,
,
,
,0,
1720 43 52632 2322 84065 26703
3
875 23 264
125 43
y, ,0,
, ,
,0,
,
40 2244 72 1955 11592 616
13
2375 5 12 3
yˆ ,
,0,
, , , ,0,0
160 5984 16 85 44
RKF7(8)




Fehlberg (1968)
7(8)次12回10点
注意:k12=k1
難点
評価回数が多い
 比較子の方が高次
 f=f(t)のとき、比較
子が不適格

2 2
,
27 27
1 1 1
, ,
9 36 12
1 1
1
, , 0,
6 24 8
5 5
25 25
, ,0,
,
12 12
16 16
1 1
1 1
, ,0, 0, ,
2 20
4 5
5 25
125 65 125
,
,0,0,
,
,
6 108
108 27 54
1 31
61 2 13
,
,0,0,0,
, ,
6 300
225 9 900
2
53 704 107 67
,2,0,0,
,
,
, ,3
3
6 45
9 90
1 91
23 976 311 19 17 1
,
,0,0,
,
,
,
, ,
3 108
108 135 54 60 6 12
2383
341 4496 301 2133 45 45 18
1,
,0,0,
,
,
,
, ,
,
4100
164 1025 82 4100 82 164 41
3
6 3 3 3 6
0,
,0,0,0,0, ,
, , , ,0
205
41 205 41 41 41
1777
341 4496 289 2193 51 33 12
1,
,0,0,
,
,
,
, ,
, , 0,1
4100
164 1025 82 4100 82 164 41
41
34 9 9 9
9 41
y,
,0,0,0,0,
, , ,
,
,
, 0, 0
840
105 35 35 280 280 840
34 9 9 9
9
41 41
yˆ , 0,0,0,0,0,
, , ,
,
,0,
,
105 35 35 280 280 840 840
外挿法

一種のメタ積分法:非常に強力
刻み幅に関するリチャードソンの外挿(h→0)
 単段法、等間隔内部点、可変次数
2
 誤差がh 展開される積分法→速い収束
 手間をかければ超高精度の結果が得られる



定番:調和点列のグラッグ法(並列化可能)
特殊な2階の場合:2倍の高速化
外挿法(2)

(可変次数)アルゴリズムの概要
0.
 1.
 2.
 3.
 4.
 5.
 6.

使用する積分法(=試験積分法)を決める
tの区間[0,H]を整数分割: hn=H/n
t=0からt=Hまで、刻み幅をhnとして積分
t=Hにおける結果をY(hn) とおく
得られたY(hk)群からY(0)(=解)を外挿
収束判定
nを増やして1.に戻る
外挿法の簡単な例

試験積分法が前進オイラー法の場合

終点t=Hでの積分誤差の振舞い
y  H ; h  y(H ;0)  hy1 (H )  h y2 (H ) 
2
必要なのはy(H;0)だけ
 多項式外挿:エイトキン・ネビルのアルゴリズム



補間表の最右下端が各回の最良近似解
調和点列 hn  H / n
グラッグ法

試験積分法:修正中点則
例外1(最初):前進オイラーで代用
 例外2(最終):後退オイラーとの平均値



不安定を押さえる工夫→分割数が偶数の必要性
積分誤差はh2展開
y  H ; h  y(H ;0)  h y2 (H )  h y4 (H ) 
2

4
逐次外挿法:エイトキン・ネビル

点列:倍調和点列 hn  H /  2n
グラッグ法(2)

グラッグ法の特徴
長所1:
 長所2:
 長所3:
 長所4:


可変次数(偶数)→精度向上が容易
プログラムが簡単(任意精度で可能)
1段法と等価→刻み幅制御が容易
小規模並列化で加速
折りたたみ法:例)4PUで、ほぼ4倍速
短所1: 低精度用途ではコスト高
 短所2: 頻繁に結果を出力するときもコスト高

グラッグ法(3)

一般の場合のアルゴリズム:調和点列を採用
f 0 :=f  y 0 ,t 0 
do(n=2,N,2){
h:=H/n; h 2 =h*2.0; y1:=y 0 +h*f 0 ; t1:=t 0 +h; f1:=f  y1 ,t1 
do  k=2,n  y k :=y k-2 +h 2 *f k-1; t k =t 0 +k*h; f k =f  y k ,t k 
yˆ n :=y n-1 +h*f n ; y n :=  y n +yˆ n  *0.5
外挿&収束判定}

添字を偶奇に分割→中間変数の数が大幅削減
グラッグ法(4)


特殊な2階常微分方程式の場合
試験積分法:中点則
d2 x
 f  x, t 
2
dt
例外1(最初):前進オイラー
 例外2(最終):後退オイラー
2
 積分誤差はh 展開


逐次外挿法:再びエイトキン・ネビル

点列:純調和点列 hn  H / n
グラッグ法(5)

特殊な2階の場合のアルゴリズム
f 0 :=f  x 0 ,t 0 
do(n=1,N){
h:=H/n; h1/2 :=h*0.5;
v1/2 :=v 0 +h1/2 *f 0 ; x1:=x 0 +h*v1/2 ; t1:=t 0 +h; f1:=f  x1 ,t1 
do  k=1,n-1 {
v k+1/2 :=v k-1/2 +h*f k ; x k+1:=x k +h*v k+1/2 ; t k =t 0 +k*h; f k =f  x k ,t k }
v n :=v n-1/2 +h1/2 *f n
外挿&収束判定}
グラッグ法(6)

折りたたみ計算:並列化による加速


Ito and Fukushima (1997, AJ)
発想:試験積分は並列化可能
h=H/nのときは、h=Hのn倍だけ時間がかかる
 PU毎の負荷を均等化+途中でも判定可能
 調和点列の場合は実装が容易

2PUの例: [(1+4),(2+3)],[(5+8),(6+7)],…
 4PUの例: [(1+8),(2+7),(3+6),(4+5)],…

グラッグ法(7)


グラッグ法の刻み幅制御
方針:外挿回数を固定化


局所誤差=最終外挿結果と、その前の差


3回外挿=7次、7回外挿=15次、…
局所誤差の次数例:6次(3回と2回目の差)
刻み幅制御のアルゴリズム(既述)

問:外挿7回の場合のアルゴリズムを書け
GBS法

Burlisch and Stoer (1965)


GBS=Gragg, Burlisch, Stoerの3人の頭文字
グラッグ法との違いは、ごくわずか
逐次外挿法:有理式(Burlisch & Stoer)
 点列:Burlisch列


しかし、グラッグ法+調和点列が最良

Deuflhard (1980)
エルミート積分法

エルミート補間に基づく積分公式


1段法陽公式=テイラー級数法
もしfの時間微分の計算が簡単なら有効

例1:部分的定数係数線形微分方程式
k
n

df
dn f
d
b
n
n k
 ay  b  t   n  a y   a  k 
dt
dt
k 0
 dt 
例2:線形微分方程式
 例3:古典重力多体問題(3階までなら)

エルミート積分法(2)

1段陰公式


導出法:exp(x)のパデ近似(問:導け)
2階4次陰公式
h
h2   df 
 df  
yn 1  yn   f n 1  f n         
2
12   dt n 1  dt n 

3階6次陰公式
3
2
2





h
h 2  df 
d
f
h
d
f
d
f  
 
yn 1  yn   f n 1  f n         
 2    2  
2
10  dt n 1  dt  n  120  dt n 1  dt  n 
エルミート積分法(3)

古典重力多体問題の場合

加速度の時間微分の計算コストは比較的小
 J
aK    3
J  K  rJK

 J
da K
  3
 rJK 
dt
J  K  rJK

  J  rJK v JK  

 rJK
 v JK  3  
5
rJK
J K 


J  GM J , rJK  xJ  xK , rJK  rJK , uJK  vJ  vK

Makino and Aarseth (1992)
h
h2
h3
x n 1  x n   v n 1  v n    an 1  an  
 jn1  jn 
2
10
120
h
h2
v n 1  v n   an 1  an    jn 1  jn 
2
12
da
j
dt
エルミート積分法(4)

陰公式の解法:予測子修正子法(後述)


代表的なモード:PEC, PECEC
予測子の例

1. テイラー級数法による1段陽公式
2
3
2
h
h
h
P
P
x n 1  x n  hv n  an  jn , v n 1  v n  han  jn
2
6
2

2. 2段陽公式(後述)

可変刻み幅のときは要注意
13.常微分方程式の
数値解法 II.多段法



多段(multistep)法
オイラー法の拡張の第三の方向
思想:過去の結果を再利用する
長所1:低コスト
 長所2:高次公式が構築しやすい
 短所1:高次ほど不安定
 短所2:出発値を別途準備する必要あり
 短所3:刻み幅の変更が困難

多段法

補間公式積分型多段法
一般1階:定番=アダムス法
 特殊2階:定番=カウエル法



一般多段法、対称多段法
陽公式と陰公式


陰公式の解法:予測子修正子法と陰陽法
出発値表の計算:外挿法がお勧め
補間公式積分型多段法

1ステップ数値積分+被積分関数の補間
t
dy
yn 1  yn   f  y, t  dt
 f  y, t 
dt
t
補間公式による分類
n1
n

等間隔:ラグランジュ補間→アダムス法
 不等間隔:ガウス・ルジャンドル型補間
→ガウス・ルジャンドル積分法、ロバット積分法

アダムス(Adams)法

発想:オイラー法で過去のfの情報も使う


J.C. Adams:海王星の予言で有名
p段陽公式(p次)
p 1
yn 1  yn  h  pk f n k
k 0
実用上、最も重要
 別名:アダムス・バッシュフォース(Bashforth)


p-1段陰公式(p次)
p 2
yn1  yn  h  f n 1k
k 0
予測子修正子法で活用
 別名:アダムス・ムールトン(Moulton)

*
pk
アダムス法(2)

公式の2大表現法

直接表現
p 1
yn 1  yn  h  pk f n k
k 0
p 2
過去のf値を保持
*
yn1  yn  h  pk
f n 1k
 係数が次数ごとに異なる
k 0
 丸め誤差に弱い


階差表現(お勧め)
p 1
yn1  yn  h  k k f n
fの階差を保持
k 0
 次数によらず統一的表現
p 2
* k
y

y

h

 丸め誤差に強い
 k  fn1
n 1
n

k 0
アダムス法(3)

直接表現の例






陽公式 Adams-Bashforth (AB)
2次
3次 y
4次 y
5次 y
6次 y
yn 1  yn 
h
 3 f n  f n1 
2
n 1  yn 
h
 23 f n  16 f n1  5 f n2 
12
n 1  yn 
h
 55 f n  59 f n1  37 f n2  9 f n3 
24
n 1  yn 
n 1  yn 
h
1901 f n  2774 f n1  2616 f n2  1274 f n3  251 f n4 
720
h
 4277 f n  7923 f n1  9982 f n2  7298 f n3  2877 f n4  475 f n5 
1440
アダムス法(4)

直接表現の例






陰公式 Adams-Moulton (AM)
h
 f n1  f n 
2
2次=台形則
h
y

y

5 f  8 f  f 
3次
12
h
y

y

 9 f  19 f  5 f  f 
4次
24
h
5次 y  y  720
 251 f  646 f  264 f  106 f  19 f 
6次 y  y  h  475 f  1427 f  798 f  482 f  173 f
yn 1  yn 
n 1
n
n 1
n 1
n
n 1
n 1
n
n 1
n
n
n 1
1440
n 1
n
n 1
n 1
n2
n 1
n
n
n 1
n2
n2
n 3
n 3
 27 f n  4 
アダムス法(5)

階差係数の母(generating)関数

階差表現の演算子解釈(陽公式)


k 0
k 0
yn1  yn  h  k  k f n  S  1  h  k  k D
S 1

 k 


hD  1    log 1   
k 0

k

m 1
係数漸化式
k
(問:導け)  0  1,  m  1   m  1  k
k 0
アダムス法(6)

階差係数の母関数(2)

階差表現の演算子解釈(陰公式)


k 0
k 0
yn1  yn  h  k*k f n1  S  1  h  k*k SD
S 1

  


hSD  log 1   
k 0

*
k

k
mk *
係数漸化式
m 1
1  k

*
*
(問:導け)  0  1,  m  
k 0 m  1  k
アダムス法(7)


階差係数
陽公式


全て正
陰公式

初項以外
全て負
1 5 3 251 95 19087 5257 1070017 25713
, ,
,
,
,
,
,
,
2 12 8 720 288 60480 17280 3628800 89600
26842253 4777223 703604254357 106364763817
,
,
,
,
95800320 17418240 2615348736000 402361344000
1166309819657 25221445 8092989203533249
,
,
,
4483454976000 98402304 32011868528640000
 m  1, ,
1 1 1  19 3 863 275 33953 8183
, , ,
,
,
,
,
,
,
2 12 24 720 160 60480 24192 3628800 1036800
3250433 4671 13695779093 2224234463
,
,
,
,
479001600 788480 2615348736000 475517952000
132282840127 2639651053 111956703448001
,
,
,
31384184832000 689762304000 32011868528640000
 m*  1,
アダムス法(8)

理論誤差=階差表現で打ち切った次の項


例:単振動
実用的誤差

h 
p

 yn   p  2sin
 ~  p  h 
2 

p
増分子=同次数の予測子と修正子の差
p 2


*
k
p 1
 yn  h  f n1    k   k 1   f nk   p1 f n1 p 
k 0


注1:陽公式、陰公式を問わず計算可能
 注2:修正子より増分子を計算する方が実用的

1階用多段法

一般q段公式:刻み幅一定を仮定
陽公式: j=0
 陰公式: j=1


dy
 f  y, t 
dt
q  j 1
q

k 0
k
yn 1k  h  k f n  j k
k 0
二つの特性多項式と安定性多項式Q

第1特性多項式、第2特性多項式
q
  z    k z q k
k 0
  z 
q  j 1
z
k 0
Q  z;      z     z 
k
q  j 1 k
差分方程式
q

斉次定数係数差分方程式
前提:等間隔刻み幅
 一般解=基本解の線形結合


基本解

k 0
k
yn  k  0
q
yn m   Ck km
k 0
q
  z    k z q k
特性根=第1特性多項式の根
k 0
 単根の場合:特性根の累乗
 km  zkm
 n重根の場合:
 km  zkm , mzkm , m  m 1 zkm ,

1階用多段法(2)

ゼロ安定性(zero stability)
常に f=0 でも安定して解けるか?
 不安定な例 yn1  2 yn  3 yn1  h  3 f n  f n1 



問: f=0でも、勝手に発振することを確かめよ
安定性条件 z0  1, zk  1
主根を除く全ての特性根の絶対値が1以下
 強安定:1未満、弱安定:=1


最安定:特性多項式=(z-1)zn…アダムス法
1階用多段法(3)

線形(linear)安定性

線形問題 f=y のとき、安定して解けるか?
線形安定性条件   h



Q(z;)の主根以外の根の絶対値が1以下

強安定:1未満、弱安定:=1
線形安定領域の二大決定法

根追跡法と境界追跡法
1階用多段法(4)

根位置(root locus)追跡法   h


個々のに対して、代数方程式 Q(z;)=0 の全
複素根を決定→非常に時間がかかる
境界(boundary)追跡法

安定境界関数


意味:Q=0のときのの値
   
  ei 
 e
i

複素面上で()をプロットし、閉領域(複数)の
内側か外側かを判定→簡単かつ高速
周期性関数


周期性(periodicity)関数 g    i  
周期安定性境界Hmaxの決定手法
Fukushima (1998):2階対称多段法で創始
 Evans & Tremaine (1999):1階対称に応用
 Yamamoto (2003):1階一般に拡張



境界追跡法の応用
条件:g()の虚部=0, 実部が絶対値最小

Hmax = g()の実部の最小絶対値
アダムス法(9)

第1特性多項式(p次)   z    z 1 z



p1
ゼロ安定性:最も強い(確認済み)
線形安定領域:別図を参照
周期安定境界(陽公式)
AB2=0.0, AB4~0.430, AB6~0.114
 AB8~0.0295, AB12~0.00759



陰公式は陽公式より、かなり安定
高次になるほど安定境界は小さくなる
アダムス法(10)

陽公式の線形安定領域: 2~5次
複素上半面 1
 高次ほど 0.8
不安定
0.6
 注意:
0.4
4,5次では
0.2
原点近くの
0
閉領域
-1.2

AB2
AB3
AB4
AB5
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
アダムス法(11)

陰公式の線形安定領域: 3~6次
2次は
台形則
 陽公式
より安定
 RKより
不安定

3.5
3
AM3
2.5
2
1.5
AM4
1
AM5
0.5
0
AM6
-6
-5
-4
-3
-2
-1
0
1階用一般多段法

1階用一般多段法の階差表現

p段p次陽公式
p 1
p 1
yn1  yn   ak  yn  h bk  f n
k
k 1

k
k 0
p段p+1次陰公式
p 1
p 1
yn1  yn   ak  yn  h b  f n 1
k
k 1
k 0
*
k
k
1階用一般多段法(2)

母関数の導入

注:a(z)はp-1次多項式
p 1


a  z    ak z , b  z    bk z , b  z    b z
k
k 1

k
*
k 0
k 0
母関数同士の関係(問:示せ)
z  1  z  a  z 
b z 
 1  z  log 1  z 
b z
b  z 
1 z
*
* k
k
1階用一般多段法(3)

特性根と階差係数aの関係(問:示せ)
k
p 1

 1 
p 1
  z    z  1   z  z j   z  z  1   ak 1   
 z  
k 1
j 1

p 1
p 1
1    z j   1   
j 1

p 1
p 1
   ak 1   
p 1k
k
k 1
問:特性根から階差係数を求める漸化式を導け
1階用一般多段法(4)

一般多段法公式の設計手順
1. 次数および陰陽の別を設定する
 2. 0以外の適当な数の特性根ziを選ぶ

注1: ziは実根もしくは複素共役根のペア
 注2: ziは単位円内から選ぶ(0安定性の確保)

3. 対応するyの階差係数akを求める
 4. fの階差係数bk (もしくは b*k)を求める
 5. 誤差定数Cを求める

1階用対称多段法

対称性の定義(1階微分方程式の場合)
 一般係数表現   
k
p 1k , k   p 1k



特性多項式
  z    z 2 p 1   z 1  ,
  z   z 2 p 1  z 1 
注意:自由度が大きい → 多数の公式群
対称多段法の光と影
長所:最大次数定理、時間対称系に適合
 短所:弱安定、奇妙な現象(刻み幅共鳴)

1階用対称多段法(2)

対称多段法の一例

4次4段対称陽公式(最適型) (-1<a<1)
4  8a
 8  2a

yn1  2a  yn  yn 2   yn 3  h 
f

f

f
 n n2 
n 1 
3
3



4次3段対称陰公式(標準型) (-1<b<1)
17  11b
 5  b

yn1  1  2b  yn  yn1   yn2  h 
f

f

f

f
 n1 n2 
 n n1 
12
 12

1階用対称多段法(3)

対称公式の2大形式:qは公式の自由度

最適型:偶数の特性根 1, exp(ij)
q
  z    z 2  1   z 2  2c j z  1
j 1

c j  cos j
標準型:奇数の特性根 1, exp(ij)
q
  z    z  1   z  2c j z  1
2
j 1
1階用対称多段法(4)

ダールキスト(Dahlquist)の最大次数定理


同一段数で最大次数の多段法=対称多段法
次数公式(q段)
陽公式最適型: p=2q
 陽公式標準型: p=2q-2
 陰公式最適型: p=2q+2
 陰公式標準型: p=2q
 アダムス法では次数は約半分(問:確認せよ)

予測子修正子法

陰公式の代表的解法

予測子 Predictor (P): yn+1を予測すること

定番:同種(=同じ特性多項式)の同次陽公式
修正子 Corrector (C): 陰公式そのもの
 評価子 Evaluator (E): fn+1を計算すること


代表的な組み合わせ(mode): PECE
P → E → C → E → 次のステップ
 陽公式の倍のコスト:どちらが良いかは難解

予測子修正子法(2)

予測子修正子法の2大モード列
PEC列:P(EC)n 、安定性に欠ける
n
 PECE列:P(EC) E 、安定性が良好


普通は収束判定せず反復回数を固定


注意:安定性は陰公式より少し悪くなる
収束するまで繰り返すのはコスト高

例外:「硬い」方程式の場合など
予測子修正子法(3)

安定性多項式:詳細は三井(1985)
QPEC  z;    z p    z     z       z   *  z     z   *  z  
PEC
 PECE


QPECE  z;      z     z    0    z     z  
安定性関数(Q=0の解)

  z  p  *   * 

PEC  
PEC

PECE PECE 
2
  0  
   0  
2
 4 0
予測子修正子法(4)

増分子(increaser):(修正子)ー(予測子)


予測子
修正子
(P)
n 1
y
 yn   ynk  h 
k 1
増分子
f
k 0
yn(C)1  yn(P)1  k(I) ynk  h k(I) f n1k
k 1

k 0
(P)
k
n k
(C)
(C)
yn(C)

y


y

h

 k nk  k fn1k
1
n
k 1

(P)
k
k 0
(I)
(C)
(P)
(I)
(C)
(I)
(C)
(P)





,



,





増分係数 k
k
k
0
0
k
k
k 1
予測子修正子法(5)

増分子が有効←右辺のy-依存性が軽い時


軽い依存性


右辺の分割(例:線型方程式、変分方程式)
dy
dy
 A t  y  b t 
 f 0  y, t   f1  t 
dt
dt
f0(y,t)の計算コストがf1(t)のコストより小さい
好例:(地球の章動など)強制振動
陰陽法

依存が非対称な場合の陰公式の活用法


Yamamoto (2003)
非対称な常微分方程式系
dw
dz
 f  w, t  ,  g  z , w, t 
dt
dt

良い例:一般2階微分方程式
2
d x
dv
 dx  dx
 a  x, , t  
 v,
 a  x, v, t 
2
dt
dt
dt
 dt 
陰陽法(2)

発想:陽公式と陰公式の組み合わせ
他に依存する部分を、先に陽公式で
 他に依存しない部分を、次に陰公式で


2階常微分方程式、アダムス法の場合
p 1
vn1  vn  h  k k an
k 0

p 1
xn1  xn  h  k*k vn1
k 0
特徴:全体では陽公式、安定性が強化

並列化するときは要注意
エルミート積分法(5)

アダムス型エルミート積分公式の導出法
1. fのエルミート fn  s    f nk Pk  s     df  Qk  s 
k 0
k 0  dt  n  k
補間公式を導出
1
yn1  yn   f n  s  ds
 2. 補間公式を積分
0


例:2階2段4次陽公式(一定刻み幅を仮定)
h
h2
yn 1  yn    f n  3 f n1   17 f n  7 f n1 
2
12

用途:既出の2階4次陰公式に対する予測子
特殊2階用多段法


d2 x
 f  x, t 
2
dt
1階用と同じ記号を使うので注意
q
q  j 1
一般q段公式
2
 k xn1k  h  k fn j k


陽: j=0、陰: j=1
k 0
k 0
二つの特性多項式と拡張特性多項式
q
  z    k z
q k
  z 
q  j 1
z
k 0
k 0
Q  z;      z      z 
2
k
q  j 1 k
補間公式積分法(2)

1ステップ数値積分+被積分関数の補間
t
dy
yn 1  yn   f  y, t  dt
 f  y, t 
dt
t
補間公式による分類
n1
n

等間隔:ラグランジュ補間→アダムス法
 不等間隔:ガウス・ルジャンドル型補間
→ガウス・ルジャンドル積分法、ロバット積分法

カウエル(Cowell)法

特殊2階常微分方程式専用
Störmer (1907) オーロラ中の荷電粒子の運動
 Cowell & Crommelin (1910) ハレー彗星の軌道



発想:2階微分を中心差分の2乗で近似
陽公式の例


2
x

2
x

x

h
fn
3点2次対称 n1
n
n1
陰公式の例

h2
3点4次対称 xn1  2 xn  xn1   f n1  10 f n1  f n1 
12
カウエル法(2)

陽公式:シュテールマー(Störmer)

一般表現と階差表現
p 1
xn1  2 xn  xn1  h2   pk f nk xn1  2 xn  xn 1  h 2
k 0

p 1
k


 k fn
k 0
陰公式:(狭義の)カウエル
p 1
p 1
k 0
k 0
2
* k
*
xn1  2 xn  xn1  h2   pk
f n1k xn1  2 xn  xn1  h  k  f n 1
カウエル法(3)

階差係数の母関数:陽公式
1
2
S

2

S

k





k
2
2 2
h D
k 0
1    log 1   


係数漸化式(問:導け)
hm1k k
 0  1,  m  1  2
k 0 m  2  k
m 1
m
1
hm  
k 1 k
カウエル法(4)

階差係数の母関数:陰公式
1


S

2

S

* k
 k 



2
2
h SD
k 0
 log 1    


係数漸化式(問:導け)
*
h

 0*  1,  m*  2 m1k k
k 0 m  2  k
m1
2
カウエル法(5)

階差係数
陽公式
1 1 19 3 863 275 33953 8183
, ,
, ,
,
,
,
,
12 12 240 40 12096 4032 518400 129600

3250433 4671 13695779093 2224234463
,
,
,
,
53222400 78848 237758976000 39626496000
 全て正
132282840127 2639651053 111956703448001
(第2項以外) 2414168064000 , 49268736000 , 2134124568576000 ,

陰公式
全て負
(第1,3,4項
以外)

 m  1,0,
1
1 1 221 19 9829 407
,0,
,
,
,
,
,
,
12 240 240 60480 6048 3628800 172800
330157 24377
4281164477
70074463
,
,
,
,
159667200 13305600 2615348736000 47551795200
1197622087 97997951
264713507083
,
,
,
896690995200 80472268800 237124952064000
 m*  1,  1,
カウエル法(6)

より良い実装法:1階階差の導入

陽公式
xn 1  xn
un 1 
h
p 1
un1  un  h k  k f n , xn 1  xn  hun 1
k 0

陰公式
p 1
un1  un  h  f n1 , xn1  xn  hun1
k 0

*
k
k
問:1階階差の導入の利点を二つ上げよ
カウエル法(7)

速度公式:1階階差を利用
エネルギーの計算などに必要
 陽公式と陰公式

p 1
p 1
vn1  un1  h k  f n
vn1  un1  h k*k f n1
k
k 0
k 0

階差係数の母関数(問:導け)

 z
k 0
k
k

 z  log 1  z 
1  z  log 1  z  

2
 z
k 0
* k
k

 z  log 1  z 
log 1  z  
2
特殊2階用多段法(2)

ゼロ安定性
主根(=重根1)を除いて、全ての特性根の絶
対値が1以下
 強安定:1未満、弱安定:=1


最も安定

特性多項式=(z-1)2zn…カウエル法
14.常微分方程式の
数値解法 III.その他







軌道運動の誤差成長
シンプレクティック(symplectic)積分法
多様体補正(manifold correction)
変数変換の重要性
変分方程式
非差分解法(ピカール法、線形変分法)
境界値問題
数値積分の誤差測定


解析解との差:非現実的
往復誤差:コスト高、非可逆系に適用不能


高次公式との差:数値不安定に弱い



終期値から逆算して初期値に戻るか?
埋込型解法(DOPRI,外挿法など)は重宝
高精度解法との差:安易、コスト高
本格的:刻み幅を半分にして比較
軌道運動の誤差成長

(摂動)二体問題を直交座標で数値積分


問:以下の経験法則を確かめよ
打切り誤差の成長(特殊な場合を除いて)
位置(特に軌道経度方向)誤差は時間の2次
 エネルギーなど他の成分は時間の1次


丸め誤差の成長 (Brouwer 1937, AJ)
位置(特に軌道経度方向)誤差は時間の1.5次
 エネルギーなど他の成分は時間の平方根

シンプレクティック積分法


シンプレクティック積分法=正準変換積分法
 H 
仮定1:ハミルトン系
dq H dp

, 



=運動方程式が正準形 dt
p dt


 q 
仮定2:ハミルトニアンが分離可能
H  T  p U  q
 多くの力学系で実現
特徴:劇的な誤差減少
打切り誤差:位置は1次成長、その他は一定
 中心力では角運動量の積分誤差=厳密に0

シンプレクティック
積分法(2)

1次1回公式(添字の微妙な違いに注意)

座標先行:Sp(h)Sq(h)
 T 
 U 
qn1  qn  h 
 , pn1  pn  h 

 p n
 q n1

運動量先行:Sq(h)Sp(h)
 U 
 T 
pn1  pn  h 
 , qn1  qn  h 

 q n
 p n1
シンプレクティック
積分法(3)

2次1回公式=1次公式の組み合わせ

座標先行:Sq(h/2) Sp(h) Sq(h/2)
qn1/ 2

運動量先行:Sp(h/2) Sq(h) Sp(h/2)
pn1/ 2

 U 
h  T 
h  T 
 qn  
, qn1  qn1/ 2  
 , pn1  pn  h 


2  p n

q
2

p

n1/ 2

n1
 T 
h  U 
h  U 
 pn  
, pn1  pn1/ 2  
 , qn1  qn  h 


2  q n

p
2

q

n1/ 2

n1
吉田の理論 (Yoshida 1990 Phys. Rev. A)
シンプレクティック
積分法(4)

4次3回公式(Ruth):座標先行版

Sq(c1h)Sp(a1h)Sq(c2h)Sp(a2h)Sq(c2h)Sp(a1h)Sq(c1h)
a1 

1
a1
1
,
a

1

2
a
,
c

,
c

 c1
2
1 1
2
3
2
2
2 2
6次7回公式(Yoshida):座標先行版の一つ

Sq(d1h)Sp(b1h)Sq(d2h)Sp(b2h)Sq(d3h)Sp(b3h)Sq(d4h)Sp(b4h)Sq
(d4h)Sp(b3h)Sq(d3h)Sp(b2h)Sq(d2h)Sp(b1h)Sq(d1h)
b1  0.784513610477560, b2  0.235573213359357, b3  1.177679984178870,
b4  1  2  b1  b2  b3  , d1 
b b
b b
b1
b b
, d 2  1 2 , d3  2 3 , d 4  3 4
2
2
2
2
シンプレクティック
積分法(5)

混合変数シンプレクティック積分法(MVS)


仮定3:摂動ハミルトニアン



Kinoshita, Nakai, and Yoshida (1991, CMDA)
H  H0  J   H1  p, q 
ほとんど全ての力学系で実現
2通りの変数空間:要素(J,)と座標(p,q)
発想:2つの変数空間で交互に正準変換

S(0)(h/2)S(1)q(h/2)S(1)p(h)S(1)q(h/2)S(0)(h/2)
多様体補正

発想:系の保存量が一定である多様体上
に載るように数値解を補正する
(p,q)
H(p,q)=E, L(p,q)=J


Nacozy (1971), Murrison (1989)
欠点:系の保存量がないと適用不可能
多様体補正(2)

発想:準保存量一定多様体に変更

Fukushima (2003-2005)
(p(t),q(t))
H(p,q)=E(t), L(p,q)=J(t)


長所:劇的な誤差減少、積分法に非依存
コスト:準保存量の時間発展も同時に追跡
多様体補正(3)

アルゴリズム
 1. 状態変数(x, v)と準保存量を同時に積分


2. 準保存量の定義式を満たすように、状態変数
を幾何学的に厳密補正


例:エネルギーE, 全角運動量L, 角運動量Z成分LZ
補正法:スケール変換、回転、射影、…
ポイント:積分の各ステップで補正
多様体補正(4)

d
x
d
v



例:摂動二体問題
 v,    3  x  a
dt
dt
r 
 運動方程式
v    dK
va
K
 
 二体エネルギーK
2  r  dt
 単スケール変換: (x, v) → (sx, sv)
 補正量決定方程式
 v2  3

 ニュートン法
  s  Ks     0
2 
r


 初期値:s=1

準保存量
2
変数変換の重要性


変数変換→微分方程式系の数値安定化
例1:摂動二体問題の正則化(regularization)
(非線形)ケプラー問題→(線形)調和振動子
 Euler, Levi-Civita, Kustaanheimo-Stiefel, …


例2:要素の変化方程式 H  H0  J   H1  J , 
作用変数J+角変数
 無摂動解=直線運動


 J  J 0  J1  t 

  0  n0t  1  t 


 厳密な数値解が得られる
変分(variational)
方程式


常微分方程式の解y、パラメータp
変分方程式=(y/p)の常微分方程式
元の常微分方程式+初期条件を偏微分して得る
 一般には非斉次線形微分方程式


例:運動方程式から派生する最小2乗法

目的関数の例 F  p   1 w  g  t ; y  p    g 2
 m m
m
2

m
(y/p)は傾斜ベクトルなどの計算に必須
変分方程式(2)

パラメータの2大分類

解の初期値y0、モデル・パラメータ
常微分方程式 dy  f y, t;  , y  y t  t


0
0
dt
と初期条件
 y 
 変分方程式 d  y   f   y 

   
, 
   t  t 
dt  y   y   y 
y 

と初期条件

jk
0
0
0
0
jk
d  y   f   y   f   y 

   

,
  0  t  t0 
dt     y          
ピカールの逐次近似法


発想:常微分方程式の数値解より、既知関数
の数値積分のほうが簡単
近似解から出発して逐次、数値積分を実行
dy
 f  y, t  , y  y0  t  t0 
dt
y
 n 1
 t   y0   f  y  s  , s  ds
t
 n
t0
あまり数値解法として使われない
 長所1:並列化で大幅な加速
 長所2:(多点)境界値問題にも対応可能

線形変分法

発想:解を基底関数の線形和で表現



y  t    k ck k  t 
変分原理、微分方程式、初期条件・境界条
件等を係数の連立方程式に翻訳
連立方程式を解く


非線形だと面倒
コツ:方程式が簡単になる直交関数系が便利
欠点:非線形性が強いと破綻
ピカール・チェビシェフ法

ピカール法+線形変分法+チェビシェフ展開



詳細: Fukushima (1997a,b, AJ)
0. 一次変換で基本区間[-1,1]に変換
1. 解と右辺を、それぞれチェビシェフ展開
y  t    k 0 Yk Tk  t 


f  y  t  , t    k 0 Fk Tk  t 

2. 微分方程式の翻訳(主要部分)
2Yk  Fk 1  Fk 1
 問:示せ
ピカール・
チェビシェフ法(2)



3. 初期・境界条件の翻訳 Y0  y0   k 1 Yk Tk  t0 

N 
4. 打ち切り   
 n
 n
5. 逐次近似アルゴリズム
Yk  y  t 
 1) 近似解の構築
n
n
y  t   f  y  t  , t 
 2) 右辺の評価
 n
 n
 3) 右辺のチェビシェフ展開 f  y  t  , t   Fk
 4) 積分の実行
 n
 n1
Fk , y0  Yk
 5) 収束判定

n
常微分方程式の
境界値問題 dy

境界値問題とは
dt
 f  y, t 
固定境界値:例 y t0   y0 , y t1   y1
t, ay t   by t   c
 自由境界値:例


代表的な方法

射的(shooting)法:多数の初期値問題に変換

うまく行かないことが多い
緩和法:差分化→大規模連立1次方程式
 変分法:直交関数展開→係数の最適化

15. 計算のテクニック



まず「高信頼性化」(後述)
何度も使うなら「高速化」
必要に迫られたら「高精度化」
誤差の分類と誤差伝播
 丸め誤差
 エンケの方法
 倍長演算

計算の高速化

プログラム全体の高速化は不必要


プロファイリング(profiling)の重要性
高速化の手法
計算機の高速化、専用計算機の利用
 より高速なアルゴリズムの採用
 近似の活用←必要な計算精度の見直し
 低精度・高精度計算のハイブリッド化
 最後の手段:実行文の高速化

プロファイリング

99:1の法則



1%の部分が99%の時間を消費する
極意:大消費時間部分だけ集中して高速化
profiling=コードの部分毎の実行時間測定
経過時間を実行回数で代用することが多い
 動的プロファイリング:実行して計測
 静的プロファイリング:机上で解析

動的プロファイリング


既存プロファイラーの活用: gprof, gcov
動的プロファイリングの実際
0. プログラムの構造化(サブルーチン、関数)
 1. 呼び出し回数の計測
→大消費時間ルーチン候補の選定
 2. 候補ルーチンごとに実行時間を測定
 3. 消費時間順にルーチンを高速化
 4. 確認:プログラムの実行時間を再測定

呼び出し回数の計測

簡単で効果的な動的プロファイリング
0.
 1.
 2.
 3.



ルーチンの全てにカウンターを設定
全体の実行の最初にカウンターをリセット
ルーチン実行時にカウンターを加算
プログラム終了時にカウンター内容を出力
注1:複数の入り口設定は、なるべく避ける
注2:カウンター加算はルーチンの最初で
実行時間の測定


準備:時間計測プログラム or ルーチン
測定の手順:「呼び出し回数計測」に準拠
注1:一回の実行時間は短いことが多い
 注2:マルチタスクOSで正確な時間計測は疑問


測定法(他の計算負荷を最小化して実行)
0. 当該ルーチンだけ実行するプログラムを準備
6
8
 1. 十分な (10 ~10 )回数だけ実行して計測
 2. 上記を多数回繰り返し、最小値を採用

静的プロファイリング

巨大プログラムでは事実上、不可能


動的プロファイリング→対象ルーチンを選別
静的プロファイリングのコツ
安直だが確実:プログラムの達人に依頼
 目標1:実行文よりループ
 目標2:多重ループなら最内側ループ
 条件分岐→動的プロファイリングを実施
 最終段階:実行文の高速化

実行文の高速化

基礎知識:基本演算の実行速度の経験則


問:各自、体感せよ
その他の実行文の実行速度の経験則
添字計算は非常に速い
 条件分岐・ループ・ルーチン呼び出しは速い
 大規模配列の呼び出しは遅い
 入出力:2進10進変換は非常に遅い
→中間ファイルは書式なしで読み書き

基本演算の実行速度


整数演算は浮動小数点演算より格段に速い
実数の加・減・乗算は、ほぼ同じ計算速度





加減乗算の回数=FLOPS
(chipにもよるが)x+y*zの形式は速い
(賢いコンパイラでは)定数式は予め計算
除算は遅く、平方根はもっと遅い
特殊関数(sin,expなど)は、数十倍遅い
実行文の高速化(2)

計算回数の軽減:中間結果の再利用
少数回の計算:速度向上の必要性は低い
 定数計算:一度だけ実施


簡単で効果的な処方箋
複合分数式→整理して単純分数式
 低次の累乗→積の繰り返し
 多項式の計算→ホーナーの方法
 複素数演算→実数ペアの演算

実行文の高速化(3)

工夫前 f:=(5.0/2.0)*sin(x)**1.5/(1.0/x+x**3)
g:=(5.0/2.0)*sin(x)**1.5/(1.0/x-x**3)

工夫後 s:=sin(x); r:=sqrt(s)
r3:=s*r; h:=2.5*r3*x
x2=x*x; x4=x2*x2
f:=h/(1.0+x4)
g:=h/(1.0-x4)
計算の高精度化

高精度化は、本当に必要か


コスト・パフォーマンスを、よくよく考えよ
高精度化の手法
計算環境の高精度化
 数学ソフトウェア
 定数表・数学ライブラリの高精度化
 アルゴリズムの見直し
 丸め誤差の防止

演算精度の測定法

簡単で確実:より高い精度の計算と比較
好例:倍精度の結果を4倍精度の結果と比較
 問:基本関数(sinなど)の計算精度を測れ



複雑だが自己完結


計算精度は引数の値に依存することに注意
保存量、恒等式、加法定理などを検査
実用的(離散計算の場合)

刻み幅を半分にして共通格子点で比較
演算精度の向上策


整数で計算できるものは整数演算で
常に倍精度計算(不足なら4倍精度)




注意:数定数の精度明示(1.0D-7など)
1程度の量に変換:スケーリング、無次元化
なるべく小さい量だけを扱う:エンケ法
丸め誤差に強い公式・アルゴリズムを採用

二大丸め誤差:桁落ちと情報落ち
桁落ちと情報落ち

桁落ち:最も危険
値がほぼ同じ数の差
 有効桁数の著しい減少


情報落ち:計算の無駄
72865.4
 72823.5
41.9
72865.4
 絶対値の大きさに大きな

32.1583
違いがある数の加減
72897.6???
 桁そろえ→小さい量の計算精度がロス
桁落ちと情報落ち(2)

問:左辺を数値計算し、右辺と比較せよ
注:無限級数の場合は、有限和の上限をいろい
ろ変えて、収束する様子を観察せよ
 どうしたら結果がより良く一致するか、工夫せよ

1 1
1  
2 6

1
 e
n 0 n !
1
2
1  12  12 
10
10
1 1
1  
2 3
106
1


 log 2
n 0 n  1

1
n
 6  1   12  1.5000005
10
n 1 10
n
誤差の分類

測定(measurement)誤差
系統(systematic)誤差:理論、機械、個人
 偶発(spontaneous)誤差:過失、自然


計算(computational)誤差:極力、減らすべし
入出力(I/O)誤差:入力時のミス、2進10進変換
 系統誤差:定数、関数、打切り、離散化
 偶発誤差:収束判定、累積(統計的意味で)
 丸め(round-off)誤差:型変換、桁落ち、情報落ち

誤差の例

2進10進変換誤差、型変換誤差
10進有限小数 2進有限小数
-n
-n
 小さい数は10 より2
 型変換での誤差例:倍精度0.1 単精度0.1


定数誤差(光速度など物理定数も要注意)


、直交多項式(Pnなど)のゼロ点、…
関数誤差(教訓:ライブラリを過信しない事)

関数ライブラリの誤差は20~100
演算の誤差伝播

絶対誤差と相対誤差 x  x  x  x 1  x 


真の値xと計算機内での近似値<x>
乗算
x  y   x 1   x     y 1   y 
 x  y  1   x   y   
局所丸め誤差, ||< 
 相対誤差
  x  y   x   y 


除算の相対誤差  x  y   x   y  



問:示せ
演算の誤差伝播(2)

加減算
x  y
  x 1   x     y 1   y  
相対誤差(問:示せ)
x x  y y   max  x , y , x  y 
  x  y 
x y
 x,yが同符号、同程度なら問題ない



x y
max  x , y     x  y 
桁落ち
情報落ち x
y    x  y   x 
y
演算の誤差伝播(3)

他の演算:結局は四則演算に帰着
関数内の誤差=原理的に推定可能
 複雑(判定を含む)な演算→正確な推定は困難

 xf   x  
 f  x  
 x 
 f  x 

関数の誤差

相対誤差伝播の感度
xf   x 
 教訓:関数のゼロ点付近は要注意
f  x
丸め誤差






有限桁計算の宿命
丸めの方式
二大丸め誤差:桁落ちと情報落ち
累積(accumulated)丸め誤差
区間(interval)演算と区間解析:非実用的
丸め誤差の軽減策
丸めの方式

最小単位ULP (Unit in the Last Place)


値xを表現するときの最終桁の1, ULP ~ |x|
四大丸め方式: 絶対誤差 x|x|
RN: 0捨1入, |x|<0.5ULP, もっとも望ましい
 RP: 切り上げ, 0<x<ULP
 RM: 切捨て, -ULP<x<0
 RZ: 0方向丸め, |x|が小さくなるように丸める


x>0ならRM, x<0ならRP
累積丸め誤差

仮定1:局所(local)丸め誤差は互いに独立
 j  k   2 jk

仮定2:局所丸め誤差の平均はゼロ
k  0

結論:累積誤差の推定値
N

k 1

k
 N
注意:実際は仮定が成立しない場合が多い
丸め誤差の軽減策



基本:単精度計算は使わない
安直:高精度計算(4倍精度、任意長演算)
方策1: 加減算の工夫


方策2: 微小量に着目


小さい量から加える、同符号でまず加算、…
エンケの方法、差分演算子、式の変形
方策3: むやみに刻み幅を小さくしない
N
a
加減算の工夫

工夫1: なるべく絶対値をそろえる

例:2分木(B-tree)和:実現が面倒な割に効果小
 a  a    a
1

2
3

 a4      a5  a6    a7  a8   
工夫2: 絶対値が増大する方向に加算


n0
n
例:絶対値が単調減少→逆順に加算 an  an1
工夫3: まず同符号で加算+最後に減算

例:交代級数
 N
  N

 1 an    a2n     a2n1 

n 1
 n1
  n1

2N
n
加減算の工夫(2)

工夫4: 加減算だけ、より高い精度で演算


難点:コスト高
決定版: デッカー(Dekker)の倍長和
加減算だけ事実上の倍精度で実施
z1:=s1 +y
 アルゴリズムの例

総和(s+=y)で、sを2変数化 z 2 :=   s1 -z1  +y  +s 2
 上位s1,下位s2 ,積み残しz2
s1:=z1 +z 2
 括弧()の有無が重要
s 2 :=  z1 -s1  +z 2

エンケ(Encke)の方法


本来:彗星の高速軌道計算 (Encke 1854)
発想:小さい量だけを扱う
仮定:問題=解ける問題+微小変化
 「解ける問題」の解析解を活用
 微小有限変化量だけを変数→桁落ちの回避



例1: 非線型方程式の解法
例2: 常微分方程式の数値解法
エンケの方法(2)

例1:ケプラー方程式の差分解法
E  e sin E  M
仮定1: 類似解が既知 E0  e0 sin E0  M 0
 仮定2: パラメータの変化量が直接与えられる
 有限変化量xの方程式 e  e  e0 , M  M  M 0

x  e cos  E0  x  sin x  N

注:xは小さい x  E / 2, N   M  e sin E  / 2
0
エンケの方法(3)

例2:摂動調和振動子の数値積分
2
d x
  0    x, t   x    x, t 
2
dt
2
d
x
 無摂動部分の運動方程式
 0 x  0
2
dt
 有限変化量xの方程式
2
d x
 0 x
2
dt
   x  t   x, t     x  t   x, t   x  t   x 
倍長倍精度演算


元祖:デッカー、単精度で倍精度を実現
発想:倍精度演算だけで4倍精度計算を実現
変数を倍精度変数のペア(上位、下位)で表現
 通常の4倍精度より3~10倍程度高速




四則演算: Dekker (1971, Numer. Math.)
平方根: Fukushima (2001, AJ)
指数関数、三角関数: Fukushima (2003)
倍長倍精度演算(2)



倍長表現 x1 , x2
和条件と正規化条件 x  x1  x2 , x2   x1
ルーチン群(講師から入手可能)
1.
 2.
 3.
 4.
 5.

入出力:read, write
加減乗算:norm, add, sq, mul
基本初等関数:pow, exp, sin, cos
派生初等関数:tan, sinh, cosh, tanh
逆関数:div, sqrt, cbrt, log, atan2, atanh
16. 計算の高信頼性化


高信頼性(reliability)化の手法
良いプログラムの書き方
プログラムの極小化
 予防(preemptive)プログラミング
 試作(prototyping)プログラミング
 リサイクルとテンプレートの活用
 プログラムの可読化・可搬化・文書化


実践的デバッグ、その他の教訓
高信頼性化の手法


正しくなければ、意味がない
高信頼性化の直接手法
計算結果の精度の確認
 入力データ・定数の妥当性の確認
 他と比較:手計算、他の手法、他人の計算
 物理的直感との比較


高信頼性化の間接かつ強力な手法

良いプログラムを書く
プログラムの極小化


巨大プログラム=エラーの温床
個々のプログラム単位は、なるべく小さく


理想:一画面で全体が一望できるサイズ
単能プログラムを組み合わせて利用
プログラムの万能化を避けよ
 UNIX流の考え方:フィルターとパイプライン


前処理・後処理・検査・入出力などの分離
予防プログラミング

プログラミングにおけるマーフィーの法則
可能性が少しでもあれば、必ずエラーは起きる
 典型:添字や引数の範囲逸脱(out of range)
 特にループ中の境界値(=最初と最後)に注意


予防(preemptive)プログラミング
「転ばぬ先の杖」
 例1: a(i) → a(max(1,min(i,IMAX)))
 例2: sqrt(x) → sqrt(abs(x)+EPSILON)

予防のコツ

名前
紛らわしい変数名の不使用(悪例:x1とxl)
 定数、変数、ルーチン名の容易な区別化


例:定数は全大文字、変数は小文字主体
意味のある名前付け(例:PI_HALF)
 類似名で違う部分は先頭に



例:positionx, positiony → x_pos, y_pos
大文字・特殊記号の活用(例:SumX, x_mean)
予防のコツ(2)

(広い意味での)型(type)

全ての変数の型を明示的に(explicit)指定

Fortranでは、implicit none宣言を忘れずに
お勧め→名前による型の(自己流)区別
 一例:(Fortran風の)先頭文字による型区別

整数:I,J,K,N、実数:R,F、複素数:W,Z
 単精度:S、倍精度:D、4倍精度:Q
 論理:L、文字:C、文字列:A、2進数:B、ポインタ:P
 単位ベクトル:U、ベクトル:V、行列:M

予防のコツ(3)

入出力
入力の安全化:自由欄(column-free)入力
 生の(raw)入力の保存 or 直接(direct)出力
 入力変数の入念な検査



検査出力:自作ルーチン(my_print等)の活用


許容範囲、桁間違い、単位間違い、型制限
検査出力専用ファイル、十分な事前検査と適合化
正規出力の整形(formatting)

グラフ化、桁そろえ、余分な情報の抑制
予防のコツ(4)

演算前の検査と適合化(adaptation)
添字(index)、除数(divisor)、引数(argument)
 検査による警告+例外処理



既定(default)値による強制適合化


例:if(x<0.0) then call ALERT
例:g=max(gMIN,min(g,gMAX)))
微小量加算による適合化

例:z=x/(y+sign(y)*EPSILON)
予防のコツ(5)

関数とサブルーチン

典型的エラー:引数の数、順序や型の不一致


引数の3大区別(注釈化するのが賢明)


効果的:(grepによる)定義・呼び出し一覧表の作成
1) 入力(=参照)、2) 出力、3) 入出力兼用
お勧め:実行状態を関数値として返す関数
実行状態:正常終了、異常終了(その1、その2…)
 本体の構成が簡易化&例外処理に便利

試作プログラミング

設計時の試作
アルゴリズムを(自己流)疑似言語で記述
 Mathematicaなどで試作アルゴリズムを試験


コーディング時の試作はトップ・ダウンで
なるべく既成品を活用←道具箱
 難しい部分は後回し=埋め草(stab)化


試作(Proto),実用(Practice),可搬(Portable)

部品を一つずつ試作・検査→固定
ルーチンのリサイクル


最も有効な生産性の向上策
名人への近道=道具箱(toolbox)の整備
自分で経験した優秀ルーチンを収集・整理
 定評あるライブラリから拝借&自己変形
 理想は自家製ライブラリ化 (mylib.h)



数定数の名前付き定数化←再利用可能
実用的な知恵:ラッピングとマクロ
ルーチンのリサイクル(2)

ラッピング(wrapping):別名「薄皮饅頭」
=既存ルーチンに、その時の用途に合うよう
に前処理・後処理をかぶせること
 引数の変換→既存ルーチン(群)の呼び出し
→結果を逆変換 (場合分けも含まれる)


道具(mini-tool)としてのマクロ
1行程度のミニルーチン→文字置換が効果的
 使いすぎに注意

テンプレートの活用

プログラム・テンプレート(template)
理想:名人の作った簡単なプログラム
 次善:定番プログラムを自分でテンプレート化
 実際:以前に作ったプログラム


テンプレートの活用法


枠組みだけ利用、部品は道具箱から選択
実用的な例:7段階テンプレート

前準備・後始末、入・出力、前・後処理、本体
7段階テンプレート







1. 前準備:型・配列の宣言、ファイルの開
(open)、全変数値のリセット、定数の計算
2. 入力:入力変数の検査・適合化・直出力
3. 前処理:変数変換、係数・準定数の計算
4. 本体:計算本体
5. 後処理:逆変換、最終結果への整形
6. 出力:中間結果の保存、最終結果表示
7. 後始末:ファイルの閉(close)、例外処理
プログラムの可読化

適切かつ多すぎないコメント



意味のある変数名→予防プログラミング
構造表現のための字下げ・サブルーチン化


参考文献、変数の意味、アルゴリズムの核心
注意:「完全構造化」に、こだわらない
外部とのインターフェースの明確化

呼び出しルーチン、入力・出力変数の明示
プログラムの可搬化

可搬性(portability)の確保
プラットフォーム、OS、プログラミング言語…
 方言(dialect)・特殊ライブラリの使用を避ける


輸出(export)のための望ましい条件
0.
 1.
 2.
 3.

プログラムの可読化
試運転プログラム・使用例・入出力例の添付
作者名・問い合わせ先・メールアドレスの明記
利用手引き(manual)の添付
プログラムの文書化

文書化(documentation)は他人の為ならず


背景説明メモの重要性


「明日の君は他人も同然」
動機、主要アルゴリズム、試作品、開発日誌
簡単で効果的な整理法
1. 主要ファイルの全てをアーカイブ化
 2. 複数の写しを別々の場所・媒体に保存
 3. ファイル名一覧+保存場所→印刷&保管

実践的デバッグ


安直:既存デバッガの活用
実践的:適切な出力文の挿入
最善:自作デバッグ出力ルーチンの呼び出し
 出力:呼び出し箇所、変数名、変数値、…


効用
1. 整形された出力→バグ発見の簡易化
 2. 無駄な出力の防止→情報の精選化
 3. 終了時にコメント化→再利用が容易

その他の教訓







1に「正しさ」、2に「わかりやすさ」
コメントは中身と一致させる
自己流は避けてライブラリを活用
Know-HowよりKnow-Who
名人のプログラムを熟読玩味
分割統治(divide and rule)
入力は打ちやすく、出力は読みやすく
その他の教訓(2)







設計は一気に、コーディングは休みながら
机上デバッグは実行時デバッグに百倍優る
汎用プログラムは不用プログラム
単一言語・同一環境に固執しない
単純なプログラムは速く、かつ間違いにくい
一つの虫の発見は百匹の虫の可能性
実数の等値(equality)比較は無意味
その他の教訓(3)







ダメなプログラムは改良より書き直し
コードの高速化よりアルゴリズムの高速化
道具箱には良いルーチンだけを精選
添字が一つ違いのエラーに注意
多入力多出力より単一入力単一出力
高精度化のときは手抜き厳禁
変更の前に、まず現状保存(save & store)
その他の教訓(4)







マジック・ナンバーは名前付き定数に
(宣言文など)同一部分はinclude化
大域変数やcommonブロックは最小限に
ジャンプ命令の利用は最小限に
無理にジャンプ命令を一掃する必要はない
一つのコンパイラを過信せずに他も試す
古いプログラムでは言語定義の変遷に留意
その他の教訓(5)







ループの初めと終わりに悪魔が潜む
if文のネストはif-then-elseの繰り返しで
if-then-elseでは最後のelseに注意
読みやすいプログラムは間違いにくい
保守する立場に立ってプログラムを書く
実行文の繰り返しはルーチンに書き直す
意味のない中間変数は、なるべく使わない
参考文献




福井他; 1999、新 数値計算、共立出版
Gear; 1971, Numerical Initial Value Problems in
Ordinary Differential Equations, Prentice-Hall
Hairer, Nørsett, and Wanner; 1993, Solving
Ordinary Differential Equations I. Nonstiff
Problems, 2nd ed., Springer-Verlag
Hairer, and Wanner; 1991, Solving Ordinary
Differential Equations II. Stiff and DifferentialAlgebraic Problems, Springer-Verlag
参考文献(2)



Henrici; 1962, Discrete Variable Methods in
Ordinary Differential Equations, Wiley (邦訳:
ヘンリッチ、清水・小林訳;1973、計算機による
常微分方程式の解法 I、II、サイエンス社)
一松・戸川編;1975、数値計算における誤差、
共立出版
茨木・福島;1991、FORTRAN77最適化プログラ
ミング、岩波書店
参考文献(3)



市田・吉本;1979、スプライン関数とその応用、
教育出版
伊理・藤野;1985、数値計算の常識、共立出版
Kernighan; 1976, Software Tools, AddisonWesley. (邦訳:カーニハン、木村訳;1981、ソフ
トウェア作法、共立出版)
参考文献(4)



Kernighan and Pike; 1999, Practice of
Programming Style, Addison-Wesley. (邦訳:
カーニハン・パイク、福崎訳;2000、プログラミング
作法、アスキー)
Kernighan and Plauger; 1974, Elements of
Programming Style, Bell Telephone Lab. (邦
訳:カーニハン・プローガー、木村訳;1976、プロ
グラム書法(第2版)、共立出版)
Lambert; 1973, Computational Methods in
Ordinary Differential Equations, Wiley
参考文献(5)





Mathews and Fink; 2004, Numerical Methods:
Using MATLAB, 4th ed., Prentice-Hall
三井;1985、数値解析入門、朝倉書店
森;2002、数値解析(第2版)、共立出版
森;1984、数値解析法、朝倉書店
Nash; 1990, Compact Numerical Methods for
Computers, 2nd ed., Adam Hilger (邦訳:ナッ
シュ、有賀・加藤・宮里訳;1996、わかるコン
ピュータ数値計算、オーム社)
参考文献(6)



二宮編;2004、数値解析のつぼ、共立出版
二宮編;2006、数値解析のわざ、共立出版
Ralston and Rabinowitz; 1978, First Course in
Numerical Analysis, 2nd ed., McGraw Hill,
reprinted from Dover, 2001 (邦訳:ラルスト
ン・ラヴィノビッツ、戸田・小野訳;1986、電子計
算機のための数値解析の理論と応用、ブレイン
図書)
参考文献(7)



Stoer and Bulirsch; 1980, Introduction to
Numerical Analysis, Springer-Verlag
山内・森口・一松編;1965、電子計算機のため
の数値計算法I,II,III、培風館
山本他;1999、これだけは知っておきたい数学
ツール、共立出版
講師
福島登志夫
 自然科学研究機構
国立天文台 (NAOJ)
 181-8588 東京都三鷹市大沢2-21-1
 [email protected]
 http://chiron.mtk.nao.ac.jp/~toshio/

ダウンロード

Power Point File