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