[BACK]Return to n_wishartd-ja.texi CVS log [TXT][DIR] Up to [local] / OpenXM / src / asir-contrib / packages / doc / n_wishartd

Annotation of OpenXM/src/asir-contrib/packages/doc/n_wishartd/n_wishartd-ja.texi, Revision 1.2

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

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>