Annotation of OpenXM/src/asir-contrib/packages/doc/n_wishartd/n_wishartd-ja.texi, Revision 1.4
1.4 ! takayama 1: %comment $OpenXM: OpenXM/src/asir-contrib/packages/doc/n_wishartd/n_wishartd-ja.texi,v 1.3 2017/03/30 06:16:36 takayama Exp $
1.3 takayama 2: %comment --- おまじない ---
1.4 ! takayama 3: \input texinfo-ja
1.1 noro 4: @iftex
5: @catcode`@#=6
6: @def@fref#1{@xrefX[#1,,@code{#1},,,]}
1.4 ! takayama 7: @def@b#1{{@bf #1}}
1.1 noro 8: @catcode`@#=@other
9: @end iftex
10: @overfullrule=0pt
1.4 ! takayama 11: @documentlanguage ja
1.1 noro 12: @c -*-texinfo-*-
13: @comment %**start of header
1.3 takayama 14: @comment --- おまじない終り ---
1.1 noro 15:
1.3 takayama 16: @comment --- GNU info ファイルの名前 ---
1.1 noro 17: @setfilename asir-contrib-n_wishartd
18:
1.3 takayama 19: @comment --- タイトル ---
1.1 noro 20: @settitle n_wishartd
21:
22: @comment %**end of header
23: @comment %@setchapternewpage odd
24:
1.3 takayama 25: @comment --- おまじない ---
1.1 noro 26: @ifinfo
27: @macro fref{name}
28: @ref{\name\,,@code{\name\}}
29: @end macro
30: @end ifinfo
31:
32: @iftex
33: @comment @finalout
34: @end iftex
35:
36: @titlepage
1.3 takayama 37: @comment --- おまじない終り ---
1.1 noro 38:
1.3 takayama 39: @comment --- タイトル, バージョン, 著者名, 著作権表示 ---
1.1 noro 40: @title n_wishartd
41: @subtitle n_wishartd User's Manual
42: @subtitle Edition 1.0
43: @subtitle Aug 2016
44:
45: @author by Masayuki Noro
46: @page
47: @vskip 0pt plus 1filll
48: Copyright @copyright{} Masayuki Noro
49: 2016. All rights reserved.
50: @end titlepage
51:
1.3 takayama 52: @comment --- おまじない ---
1.1 noro 53: @synindex vr fn
1.3 takayama 54: @comment --- おまじない終り ---
1.1 noro 55:
1.3 takayama 56: @comment --- @node は GNU info, HTML 用 ---
57: @comment --- @node の引数は node-name, next, previous, up ---
1.1 noro 58: @node Top,, (dir), (dir)
59:
1.3 takayama 60: @comment --- @menu は GNU info, HTML 用 ---
61: @comment --- chapter 名を正確に並べる ---
1.1 noro 62: @menu
1.3 takayama 63: * matrix 1F1 の対角領域上への制限パッケージ n_wishartd.rr ::
1.1 noro 64: * Index::
65: @end menu
66:
1.3 takayama 67: @comment --- chapter の開始 ---
68: @comment --- 親 chapter 名を正確に ---
69: @node matrix 1F1 の対角領域上への制限パッケージ n_wishartd.rr ,,, Top
70: @chapter matrix 1F1 の対角領域上への制限パッケージ n_wishartd.rr
71: @comment --- section 名を正確に並べる ---
1.1 noro 72: @menu
1.3 takayama 73: * matrix 1F1 の対角領域上への制限::
74: * 制限した関数の数値計算::
75: * 部分分数係数の微分作用素::
76: * Runge-Kutta 法の試験的実装::
1.1 noro 77: @end menu
78:
1.3 takayama 79: このマニュアルでは, asir-contrib パッケージに収録されている,
80: matrix 1F1 が対角領域上で満たす微分方程式系を計算するパッケージ
81: @samp{n_wishartd.rr} について解説する.
82: このパッケージを使うには, まず @samp{n_wishartd.rr} をロードする.
1.1 noro 83: @example
84: [...] load("n_wishartd.rr");
85: @end example
86: @noindent
1.3 takayama 87: このパッケージの函数を呼び出すには, 全て @code{n_wishartd.} を先頭につける.
1.1 noro 88:
1.3 takayama 89: @comment --- section の開始 ---
90: @comment --- 書体指定について ---
91: @comment --- @code{} はタイプライタ体表示 ---
92: @comment --- @var{} は斜字体表示 ---
93: @comment --- @b{} はボールド表示 ---
94: @comment --- @samp{} はファイル名などの表示 ---
1.1 noro 95:
1.3 takayama 96: @node matrix 1F1 の対角領域上への制限,,, matrix 1F1 の対角領域上への制限パッケージ n_wishartd.rr
97: @section matrix 1F1 の対角領域上への制限
1.1 noro 98:
99: @menu
100: * n_wishartd.diagpf::
101: * n_wishartd.message::
102: @end menu
103:
1.3 takayama 104: @node n_wishartd.n_wishartd.diagpf,,, matrix 1F1 の対角領域上への制限
1.1 noro 105:
106: @subsection @code{n_wishartd.diagpf}
107: @findex n_wishartd.diagpf
108:
109: @table @t
110: @item n_wishartd.diagpf(@var{m},@var{blocks})
1.3 takayama 111: @var{m}変数の 1F1 が満たす方程式を, @var{blocks} で指定される
112: 対角領域上に制限した微分方程式系を計算する.
1.1 noro 113: @end table
114:
115: @table @var
116: @item return
1.3 takayama 117: @var{[E1,E2,...]} なるリスト, 各 @var{Ei} は
118: 部分分数を係数とする微分作用素で, 制限した 1F1を零化する.
1.1 noro 119:
120: @item m
1.3 takayama 121: 自然数
1.1 noro 122: @item vars
1.3 takayama 123: @var{[[s1,e1],[s2,e2],...]} なるリスト.
1.1 noro 124: @item options
1.3 takayama 125: 下の説明参照.
1.1 noro 126: @end table
127:
128: @itemize @bullet
1.3 takayama 129: @item @var{m}変数の 1F1 が満たす方程式を, @var{blocks} で指定される
130: 対角領域上に制限した微分方程式系を計算する.
131: @item @var{blocks} の各成分 @var{[s,e]} は @var{ys=y(s+1)=...=ye} を意味する.
132: このブロックを代表する変数は @var{ye} である.
133: @item @var{blocks} には全ての変数が現れるように指定する. 特に, 一つの変数からなる
134: ブロックは @var{[s,s]} と指定する.
135: @item この関数が常に成功することは証明されていないが, 現在のところ, 各変数ブロックサイズが
136: 36 以下なら成功することは証明されている.
137: @item 出力される微分作用素のフォーマットに関しては @ref{部分分数を係数とする微分作用素} を
138: 参照.
1.1 noro 139: @end itemize
140:
141: @example
142: [2649] Z=n_wishartd.diagpf(5,[[1,3],[4,5]]);
143: [
144: [[[[-1,[]]],(1)*<<0,0,0,0,3,0>>],
145: [[[-2,[[y1-y4,1]]],[-2,[[y4,1]]]],(1)*<<0,1,0,0,1,0>>],
146: [[[9/2,[[y1-y4,1]]],[-3*c+11/2,[[y4,1]]],[3,[]]],(1)*<<0,0,0,0,2,0>>],
147: ...
148: [[[-6*a,[[y1-y4,1],[y4,1]]],[(4*c-10)*a,[[y4,2]]],[-4*a,[[y4,1]]]],
149: (1)*<<0,0,0,0,0,0>>]],
150: [[[[-1,[]]],(1)*<<0,4,0,0,0,0>>],
151:
152: [[[-6,[[y1-y4,1]]],[-6*c+10,[[y1,1]]],[6,[]]],(1)*<<0,3,0,0,0,0>>],
153: [[[5,[[y1-y4,1]]],[-5,[[y1,1]]]],(1)*<<0,2,0,0,1,0>>],
154: ...
155: [[[21*a,[[y1-y4,2],[y1,1]]],[(36*c-87)*a,[[y1-y4,1],[y1,2]]],
156: [-36*a,[[y1-y4,1],[y1,1]]],[(18*c^2-84*c+96)*a,[[y1,3]]],
157: [-9*a^2+(-36*c+87)*a,[[y1,2]]],[18*a,[[y1,1]]]],(1)*<<0,0,0,0,0,0>>]]
158: ]
159: @end example
160:
1.3 takayama 161: @node n_wishartd.message,,, matrix 1F1 の対角領域上への制限
1.1 noro 162:
163: @subsection @code{n_wishartd.message}
164: @findex n_wishartd.message
165:
166: @table @t
167: @item n_wishartd.message(@var{onoff})
1.3 takayama 168: 計算中のメッセージ出力をon/off する.
1.1 noro 169: @end table
170:
171: @table @var
172: @item onoff
1.3 takayama 173: @var{onoff=1} のときメッセージを出力し, @var{onoff=0} のときしない.
1.1 noro 174: @end table
175:
176: @itemize @bullet
1.3 takayama 177: @item 計算中のメッセージ出力を on/off する. デフォルトは off である.
1.1 noro 178: @end itemize
179:
1.3 takayama 180: @node 制限した関数の計算,,, matrix 1F1 の対角領域上への制限パッケージ n_wishartd.rr
181: @section 制限した関数の計算
1.1 noro 182:
183: @menu
184: * n_wishartd.prob_by_hgm::
185: * n_wishartd.prob_by_ps::
186: * n_wishartd.ps::
187: @end menu
188:
1.3 takayama 189: @node n_wishartd.prob_by_hgm,,, 制限した関数の計算
1.1 noro 190: @subsection @code{n_wishartd.prob_by_hgm}
191: @findex n_wishartd.prob_by_hgm
192:
193: @table @t
1.2 noro 194: @item n_wishartd.prob_by_hgm(@var{m},@var{n},@var{[p1,p2,...]},@var{[s1,s2,...]},@var{t}[|@var{options}])
1.3 takayama 195: HGM により重複固有値を持つ共分散行列に対する Wishart 行列の最大固有値の
196: 分布関数の値を計算する.
1.1 noro 197: @end table
198:
199: @table @var
200: @item return
201: @item m
1.3 takayama 202: 変数の個数
1.1 noro 203: @item n
1.3 takayama 204: 自由度
1.1 noro 205: @item [p1,p2,...]
1.3 takayama 206: 重複固有値の個数のリスト
1.1 noro 207: @item [s1,s2,...]
1.3 takayama 208: 各重複固有値
1.1 noro 209: @end table
210:
211: @itemize @bullet
212: @item
1.3 takayama 213: 固有値 @var{si} を @var{pi} 個もつ対角行列を共分散行列とする Wishart 行
214: 列の最大固有値 @var{l1}の分布関数の値 @var{Pr[l1<t]} を計算する.
1.1 noro 215:
1.3 takayama 216: @item ステップ数を指定したルンゲ=クッタ法を, ステップ数を 2 倍しながら
217: 一つ前の計算結果との相対誤差が @var{eps} (デフォルトで @var{10^(-4)})
218: になるまで繰り返す.
1.1 noro 219: @item
1.3 takayama 220: @var{eq} オプション指定がない場合, @var{[p1,p2,...]} で指定される対角領
221: 域に制限した微分方程式系を計算する. 指定がある場合, オプションとして指
222: 定されたリストをチェックなしに制限した方程式と見なして計算する.
223: @item @var{eps}オプションが指定された場合, 指定された値を @var{eps} として計算する.
224: @item @var{td} オプションが指定された場合, 初期ベクトル計算のためのべき級数を @var{td} で
225: 指定された全次数まで計算する (デフォルトは100).
226: @item @var{rk} オプションが指定された場合, 指定された次数のルンゲ=クッタ法を用いる.
227: 許される値は 4 または 5, でデフォルトは 5である.
228: @item べき級数解の計算の困難さ, およびパフィアン行列の計算の困難さから, ブロック数が 2 以下の場合にのみ
229: 実用性がある.
1.1 noro 230: @end itemize
231:
232: @example
233: [...] n_wishartd.message(1)$
234: [...] P=n_wishartd.prob_by_hgm(10,100,[9,1],[1/100,1],100|eps=10^(-6));
235: ...
236: [x0=,8/25]
237: Step=10000
238: [0]
239: [8.23700622458446e-17,8.23700622459772e-17]
240: ...
241: Step=1280000
242: [0][100000][200000][300000]...[900000][1000000][1100000][1200000]
243: [0.516246820120598,0.516246820227214]
244: [log ratio=,4.84611265040128]
245:
246: Step=2560000
247: [0][100000][200000][300000]...[2200000][2300000][2400000][2500000]
248: [0.516246912003845,0.516246912217004]
249: [log ratio=,4.93705929488356]
250: [diag,18.6292,pfaffian,1.09207,ps,41.0026,rk,213.929]
251: 0.516246912217004
252: 266.4sec + gc : 8.277sec(276.8sec)
253: @end example
254:
1.3 takayama 255: @node n_wishartd.prob_by_ps,,, 制限した関数の計算
1.1 noro 256: @subsection @code{n_wishartd.prob_by_ps}
257: @findex n_wishartd.prob_by_ps
258:
259: @table @t
260: @item n_wishartd.prrob_by_ps(@var{m},@var{n},@var{[p1,p2,...]},@var{[s1,s2,...]},@var{t}[|@var{options}])
1.3 takayama 261: べき級数により重複固有値を持つ共分散行列に対する Wishart 行列の最大固有値の
262: 分布関数の値を計算する.
1.1 noro 263: @end table
264:
265: @table @var
266: @item m
1.3 takayama 267: 変数の個数
1.1 noro 268: @item n
1.3 takayama 269: 自由度
1.1 noro 270: @item [p1,p2,...]
1.3 takayama 271: 重複固有値の個数のリスト
1.1 noro 272: @item [s1,s2,...]
1.3 takayama 273: 各重複固有値
1.1 noro 274: @end table
275:
276: @itemize @bullet
277: @item
1.3 takayama 278: 直前の値との相対誤差が @var{eps} (デフォルト値は @var{10^(-4)}) 以下に
279: なるまで, べき級数を全次数ごとに計算する. その値から分布関数の値を計算
280: して返す.
281: @item @var{eps}オプションが指定された場合, 指定された値を @var{eps} として計算する.
282: @var{eq} オプション指定がない場合, @var{[p1,p2,...]} で指定される対角領
283: 域に制限した微分方程式系を計算する. 指定がある場合, オプションとして指
284: 定されたリストをチェックなしに制限した方程式と見なして計算する.
285: @item @var{t} の値が小さい場合にのみ実用的に用いることができる.
1.1 noro 286: @end itemize
287:
288: @example
289: [...] Q=n_wishartd.prob_by_ps(10,100,[9,1],[1/100,1],1/2);
290: ...
291: [I=,109,act,24.9016,actmul,0,gr,19.7852]
292: 2.69026137621748e-165
293: 61.69sec + gc : 2.06sec(64.23sec)
294: [...] R=n_wishartd.prob_by_hgm(10,100,[9,1],[1/100,1],1/2|td=50);
295: [diag,15.957,pfaffian,1.00006,ps,5.92437,rk,1.29208]
296: 2.69026135182769e-165
297: 23.07sec + gc : 1.136sec(24.25sec)
298: @end example
299:
1.3 takayama 300: @node n_wishartd.ps,,, 制限した関数の計算
1.1 noro 301: @subsection @code{n_wishartd.ps}
302: @findex n_wishartd.ps
303:
304: @table @t
305: @item n_wishartd.ps(@var{z},@var{v},@var{td})
1.3 takayama 306: 微分方程式系のべき級数解を指定された全次数まで計算する.
1.1 noro 307: @end table
308:
309: @table @var
310: @item return
1.3 takayama 311: 多項式リスト
1.1 noro 312:
313: @item z
1.3 takayama 314: 部分分数係数の微分作用素のリスト
1.1 noro 315: @item v
1.3 takayama 316: 変数リスト
1.1 noro 317: @end table
318:
319: @itemize @bullet
320: @item
1.3 takayama 321: 結果は @var{[p,pd]} なるリストで, @var{p} は @var{td} 次まで求めたべき級数解, @var{pd} は
322: @var{p} の @var{td} 次部分である.
323: @item @var{z} は, @var{v} に指定される変数以外のパラメタを含んではいけない.
1.1 noro 324: @end itemize
325:
326: @example
327: [...] Z=n_wishartd.diagpf(10,[[1,5],[6,10]])$
328: [...] Z0=subst(Z,a,(10+1)/2,c,(10+100+1)/2)$
329: [...] PS=n_wishartd.ps(Z0,[y1,y6],10)$
330: [...] PS[0];
331: 197230789502743383953639/230438384724900975787223158176000*y1^10+
332: ...
333: +(6738842542131976871672233/1011953706634779427957034268904320*y6^9
334: ...+3932525/62890602*y6^2+1025/4181*y6+55/111)*y1
335: +197230789502743383953639/230438384724900975787223158176000*y6^10
336: +...+1395815/62890602*y6^3+3175/25086*y6^2+55/111*y6+1
337: @end example
338:
1.3 takayama 339: @node 部分分数係数の微分作用素,,, matrix 1F1 の対角領域上への制限パッケージ n_wishartd.rr
340: @section 部分分数係数の微分作用素
1.1 noro 341:
342: @menu
1.3 takayama 343: * 部分分数の表現::
344: * 部分分数係数の微分作用素の表現::
345: * 部分分数係数の微分作用素の演算::
1.1 noro 346: @end menu
347:
1.3 takayama 348: @node 部分分数の表現,,, 部分分数係数の微分作用素
349: @subsection 部分分数の表現
1.1 noro 350:
1.3 takayama 351: matrix 1F1 が満たす微分方程式の係数は @var{1/yi}, @var{1/(yi-yj)} の定
352: 数倍の和として書かれている. さらに, ロピタル則を用いた対角領域への制限
353: アルゴリズムの結果も同様に部分分数の和として書ける.
1.1 noro 354:
355: @itemize @bullet
356: @item
1.3 takayama 357: 分母に現れる @var{yi0^n0(yi1-yj1)^n1(yi2-yj2)^n2...(yik-yjk)^nk} は
358: @var{[[yi0,n0],[yi1-yj1,n1],...,[yik-yjk,nk]]} なる形のリストとして表現
359: される. ここで, 各因子 @var{yi-yj} は @var{i>j} を満たし, さらに因子は
360: ある一定の順序で整列される.
1.1 noro 361: @item
1.3 takayama 362: @var{f} を上のようなべき積とし, @var{c} を定数とするとき, 単項式にあた
363: る @var{c/f} は @var{[c,f]} で表現される. @var{f=[]} の場合, 分母が 1
364: であることを意味する.
1.1 noro 365: @item
1.3 takayama 366: 最後に, @var{c1/f1+...+ck/fk} は @var{[[c1,f1],...,[ck,fk]]} と表現され
367: る. ここでも, 各項はある一定の順序で整列される.
1.1 noro 368: @item
1.3 takayama 369: 部分分数を通分して簡約した結果, 0 になることもあることに注意する.
1.1 noro 370: @end itemize
371:
1.3 takayama 372: @node 部分分数係数の微分作用素の表現,,, 部分分数係数の微分作用素
373: @subsection 部分分数係数の微分作用素の表現
1.1 noro 374:
1.3 takayama 375: 前節の部分分数を用いて, それらを係数とする微分作用素が表現できる.
376: @var{f1,...,fk} を部分分数の表現, @var{d1,...,dk} を分散表現単項式 (現
377: 在設定されている項順序で @var{d1>...>dk}) とするとき, 微分作用素
378: @var{f1*d1+...+fk*dk} が@var{[f1,d1],...[fk,dk]]}で表現される.
1.1 noro 379:
1.3 takayama 380: @node 部分分数係数の微分作用素の演算,,, 部分分数係数の微分作用素
381: @subsection 部分分数係数の微分作用素の演算
1.1 noro 382:
383: @menu
384: * n_wishartd.wsetup::
385: * n_wishartd.addpf::
386: * n_wishartd.mulcpf::
387: * n_wishartd.mulpf::
388: * n_wishartd.muldpf::
389: @end menu
390:
1.3 takayama 391: @node n_wishartd.wsetup,,, 部分分数係数の微分作用素の演算
1.1 noro 392: @subsubsection @code{n_wishartd.wsetup}
393: @findex n_wishartd.wsetup
394:
395: @table @t
396: @item n_wishartd.wsetup(@var{m})
397: @end table
398:
399: @table @var
400: @item m
1.3 takayama 401: 自然数
1.1 noro 402: @end table
403:
404: @itemize @bullet
1.3 takayama 405: @item @var{m} 変数の計算環境をセットする. 変数は @var{y0,y1,...,ym}, @var{dy0,...,dym}
406: で @var{y0, dy0} は中間結果の計算のためのダミー変数である.
1.1 noro 407: @end itemize
408:
1.3 takayama 409: @node n_wishartd.addpf,,, 部分分数係数の微分作用素の演算
1.1 noro 410: @subsubsection @code{n_wishartd.addpf}
411: @findex n_wishartd.addpf
412: @table @t
413: @item n_wishartd.addpf(@var{p1},@var{p2})
414: @end table
415:
416: @table @var
417: @item return
1.3 takayama 418: 部分分数係数の微分作用素
1.1 noro 419: @item p1, p2
1.3 takayama 420: 部分分数係数の微分作用素
1.1 noro 421: @end table
422:
423: @itemize @bullet
1.3 takayama 424: @item 微分作用素 @var{p1}, @var{p2} の和を求める.
1.1 noro 425: @end itemize
426:
1.3 takayama 427: @node n_wishartd.mulcpf,,, 部分分数係数の微分作用素の演算
1.1 noro 428: @subsubsection @code{n_wishartd.mulcpf}
429: @findex n_wishartd.mulcpf
430: @table @t
431: @item n_wishartd.mulcpf(@var{c},@var{p})
432: @end table
433:
434: @table @var
435: @item return
1.3 takayama 436: 部分分数係数の微分作用素
1.1 noro 437: @item c
1.3 takayama 438: 部分分数
1.1 noro 439: @item p
1.3 takayama 440: 部分分数係数の微分作用素
1.1 noro 441: @end table
442:
443: @itemize @bullet
1.3 takayama 444: @item 部分分数 @var{c} と微分作用素 @var{p} の積を計算する.
1.1 noro 445: @end itemize
446:
1.3 takayama 447: @node n_wishartd.mulpf,,, 部分分数係数の微分作用素の演算
1.1 noro 448: @subsubsection @code{n_wishartd.mulpf}
449: @findex n_wishartd.mulpf
450: @table @t
451: @item n_wishartd.mulpf(@var{p1},@var{p2})
452: @end table
453:
454: @table @var
455: @item return
1.3 takayama 456: 部分分数係数の微分作用素
1.1 noro 457: @item p1, p2
1.3 takayama 458: 部分分数係数の微分作用素
1.1 noro 459: @end table
460:
461: @itemize @bullet
1.3 takayama 462: @item 微分作用素 @var{p1}, @var{p2} の積を計算する.
1.1 noro 463: @end itemize
464:
1.3 takayama 465: @node n_wishartd.muldpf,,, 部分分数係数の微分作用素の演算
1.1 noro 466: @subsubsection @code{n_wishartd.muldpf}
467: @findex n_wishartd.muldpf
468: @table @t
469: @item n_wishartd.muldpf(@var{y},@var{p})
470: @end table
471:
472: @table @var
473: @item return
1.3 takayama 474: 部分分数係数の微分作用素
1.1 noro 475: @item y
1.3 takayama 476: 変数
1.1 noro 477: @item p
1.3 takayama 478: 部分分数係数の微分作用素
1.1 noro 479: @end table
480:
481: @itemize @bullet
1.3 takayama 482: @item 変数 @var{y} に対し, 微分作用素 @var{dy} と @var{p} の微分作用素としての
483: 積を計算する.
1.1 noro 484: @end itemize
485:
486: @example
487: [...] n_wishartd.wsetup(4)$
488: [...] P=n_wishartd.wishartpf(4,1);
489: [[[[1,[]]],(1)*<<0,2,0,0,0>>],[[[1/2,[[y1-y2,1]]],[1/2,[[y1-y3,1]]],
490: ...,[[[-a,[[y1,1]]]],(1)*<<0,0,0,0,0>>]]
491: [...] Q=n_wishartd.muldpf(y1,P);
492: [[[[1,[]]],(1)*<<0,3,0,0,0>>],[[[1/2,[[y1-y2,1]]],[1/2,[[y1-y3,1]]],
493: ...,[[[a,[[y1,2]]]],(1)*<<0,0,0,0,0>>]]
494: @end example
495:
1.3 takayama 496: @node Runge-Kutta 法の試験的実装,,, matrix 1F1 の対角領域上への制限パッケージ n_wishartd.rr
497: @section Runge-Kutta 法の試験的実装
1.1 noro 498:
499: @menu
500: * rk_ratmat::
501: @end menu
502:
1.3 takayama 503: @node rk_ratmat,,, Runge-Kutta 法の試験的実装
1.1 noro 504:
1.3 takayama 505: @code{n_wishartd.ps_by_hgm} では, パフィアン行列を計算したあと, 与えられたステップ数で
506: Runge-Kutta 法を実行して近似解の値を計算する組み込み関数 @code{rk_ratmat} を実行している.
507: この関数を, 値が与えられた精度で安定するまでステップ数を2倍しながら繰り返して実行する.
508: @code{rk_ratmat} 自体, ある程度汎用性があるので, ここでその使用法を解説する.
1.1 noro 509:
510: @subsection @code{rk_ratmat}
511: @findex rk_ratmat
512:
513: @table @t
514: @item rk_ratmat(@var{rk45},@var{num},@var{den},@var{x0},@var{x1},@var{s},@var{f0})
1.3 takayama 515: 有理関数係数のベクトル値一階線形常微分方程式系を Runge-Kutta 法で解く
1.1 noro 516: @end table
517:
518: @table @var
519: @item return
1.3 takayama 520: 実数のリスト
1.1 noro 521:
522: @item rk45
1.3 takayama 523: 4 または 5
1.1 noro 524: @item num
1.3 takayama 525: 定数行列の配列
1.1 noro 526: @item den
1.3 takayama 527: 多項式
1.1 noro 528: @item x0, x1
1.3 takayama 529: 実数
1.1 noro 530: @item s
1.3 takayama 531: 自然数
1.1 noro 532: @item f0
1.3 takayama 533: 実ベクトル
1.1 noro 534: @end table
535:
536: @itemize @bullet
537: @item
1.3 takayama 538: 配列 @var{num} のサイズを @var{k} とするとき,
539: @var{P(x)=1/den(num[0]+num[1]x+...+num[k-1]x^(k-1))} に対し @var{dF/dx = P(x)F}, @var{F(x0)=f0} を
540: Runge-Kutta 法で解く.
1.1 noro 541: @item
1.3 takayama 542: @var{rk45} が 4 のとき 4 次 Runge-Kutta, 5 のとき 5 次 Runge-Kutta アルゴリズムを実行する.
543: 実験的実装のため, adaptive アルゴリズムは実装されていない.
1.1 noro 544: @item
1.3 takayama 545: @var{s} はステップ数で, 刻み幅は@var{(x1-x0)/s} である.
1.1 noro 546: @item
1.3 takayama 547: @var{f0} がサイズ@var{n} のとき, @var{num} の各成分は @var{n} 次正方行列である.
1.1 noro 548: @item
1.3 takayama 549: 結果は, 長さ @var{s} の実数リスト @var{[r1,...,rs]} で, @var{ri} は @var{i} ステップ目に計算された
550: 解ベクトルの第0成分である. 次のステップに進む前に解ベクトルを @var{ri} で割るので, 最終的に
551: 解 @var{F(x1)} の第 0 成分が @var{rs*r(s-1)*...*r1} となる.
552: @item 方程式が線形なので, Runge-Kutta の各ステップも線形となることを利用し,
553: 第0成分を1に正規化することで, 途中の解の成分が倍精度浮動小数の
554: 範囲に収まることを期待している. 初期ベクトル @var{f0} の成分が倍精度浮動小数に収まらない場合
555: は, @var{f0} を正規化してから @code{rk_ratmat} を実行し, 前項の結果に @var{f0} の第 0 成分をかければ
556: よい.
1.1 noro 557: @end itemize
558:
559: @example
560: [...] F=ltov([sin(1/x),cos(1/x),sin(1/x^2),cos(1/x^2)]);
561: [ sin((1)/(x)) cos((1)/(x)) sin((1)/(x^2)) cos((1)/(x^2)) ]
562: [...] F0=map(eval,map(subst,F,x,1/10));
563: [ -0.54402111088937 -0.839071529076452 -0.506365641109759 0.862318872287684 ]
564: [...] N0=matrix(4,4,[[0,0,0,0],[0,0,0,0],[0,0,0,-2],[0,0,2,0]])$
565: [...] N1=matrix(4,4,[[0,-1,0,0],[1,0,0,0],[0,0,0,0],[0,0,0,0]])$
566: [...] N=ltov([N0,N1])$
567: [...] D=x^3$
568: [...] R=rk_ratmat(5,N,D,1/10,10,10^4,F0)$
569: [...] for(T=R,A=1;T!=[];T=cdr(T))A *=car(T)[1];
570: [...] A;
571: 0.0998334
572: [...] F1=map(eval,map(subst,F,x,10));
573: [ 0.0998334166468282 0.995004165278026 0.00999983333416666 0.999950000416665 ]
574: @end example
575:
576:
1.3 takayama 577: @comment --- おまじない ---
1.1 noro 578: @node Index,,, Top
579: @unnumbered Index
580: @printindex fn
581: @printindex cp
582: @iftex
583: @vfill @eject
584: @end iftex
585: @summarycontents
586: @contents
587: @bye
1.3 takayama 588: @comment --- おまじない終り ---
1.1 noro 589:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to n_wishartd-ja.texi CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM / src / asir-contrib / packages / doc / n_wishartd