[BACK]Return to n_wishartd-en.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-en.texi, Revision 1.4

1.4     ! takayama    1: %comment $OpenXM: OpenXM/src/asir-contrib/packages/doc/n_wishartd/n_wishartd-en.texi,v 1.3 2016/08/29 05:31:48 noro Exp $
1.1       noro        2: %comment --- おまじない ---
1.4     ! takayama    3: \input texinfo
1.1       noro        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})
1.2       noro      163: starts/stops displaying messages during computation.
1.1       noro      164: @end table
                    165:
                    166: @table @var
                    167: @item onoff
1.2       noro      168: Start displaying messages if @var{onoff}=1.
                    169: Stop displaying messages if @var{onoff}=0.
1.1       noro      170: @end table
                    171:
                    172: @itemize @bullet
1.2       noro      173: @item This function starts/stops displaying messages during computation.
                    174: The default value is set to 0.
1.1       noro      175: @end itemize
                    176:
                    177: @node Numerical comptation of restricted function,,, n_wishartd.rr
                    178: @section Numerical comptation of restricted function
                    179:
                    180: @menu
                    181: * n_wishartd.prob_by_hgm::
                    182: * n_wishartd.prob_by_ps::
                    183: * n_wishartd.ps::
                    184: @end menu
                    185:
                    186: @node n_wishartd.prob_by_hgm,,, Numerical comptation of restricted function
                    187: @subsection @code{n_wishartd.prob_by_hgm}
                    188: @findex n_wishartd.prob_by_hgm
                    189:
                    190: @table @t
1.2       noro      191: @item n_wishartd.prob_by_hgm(@var{m},@var{n},@var{[p1,p2,...]},@var{[s1,s2,...]},@var{t}[|@var{options}])
                    192: computes the value of the distribution function of the largest eigenvalue of a Wishart matrix.
1.1       noro      193: @end table
                    194:
                    195: @table @var
                    196: @item return
                    197: @item m
1.2       noro      198: The number of variables.
1.1       noro      199: @item n
1.2       noro      200: The degrees of freedom.
                    201: @item [p1,p2,...,pk]
                    202: A list of the multiplicities of repeated eigenvalues.
                    203: @item [s1,s2,...,sk]
                    204: A list of repeated eigenvalues.
1.1       noro      205: @end table
                    206:
                    207: @itemize @bullet
                    208: @item
1.2       noro      209: Let @var{l1} be the largest eigenvalue of a Wishart matrix.
                    210: Let @var{Pr[l1<t]} be the distribution function of @var{l1}.
                    211: The function @code{n_wishartd.prob_by_hgm} computes the value of the distribution function by using HGM
                    212: for a covariance matrix which has repeated eigenvalues @var{si} with multiplicity @var{pi} (@var{i=1,...,k}).
                    213:
                    214: @item This function repeats a Runge-Kutta method for the Pfaffian system by doubling the step size
                    215: until the relative error between the current result and the previous result is less than @var{eps},
                    216: The default value of @var{eps} is @var{10^(-4)}.
1.1       noro      217: @item
1.2       noro      218: If an option @var{eq} is not set,
                    219: a system of PDES satisfied by 1F1 on the diagonal region specified by @var{[p1,p2,...]} is computed.
                    220: If an option @var{eq} is set, the list specified by @var{eq} is regarded as the correct system of PDEs.
                    221: @item If an option @var{eps} is set, the value is used as @var{eps}.
                    222: @item If an option @var{td} is set, the truncated power series solution for computing the initial vector
                    223: is computed up to the total degree specified by @var{td}. The default value is 100.
                    224: @item If an option @var{rk} is set, it is regarded as the order of a Runge-Kutta method. The default vaule is 5.
                    225: @item It is recommended to use this function only when @var{k<=2} where @var{k} is the number of diagonal blocks because of
                    226: the difficulty of the truncated power series solution and the difficulty of computation of the Pfaffian matrices.
1.1       noro      227: @end itemize
                    228:
                    229: @example
                    230: [...] n_wishartd.message(1)$
                    231: [...] P=n_wishartd.prob_by_hgm(10,100,[9,1],[1/100,1],100|eps=10^(-6));
                    232: ...
                    233: [x0=,8/25]
                    234: Step=10000
                    235: [0]
                    236: [8.23700622458446e-17,8.23700622459772e-17]
                    237: ...
                    238: Step=1280000
                    239: [0][100000][200000][300000]...[900000][1000000][1100000][1200000]
                    240: [0.516246820120598,0.516246820227214]
                    241: [log ratio=,4.84611265040128]
                    242:
                    243: Step=2560000
                    244: [0][100000][200000][300000]...[2200000][2300000][2400000][2500000]
                    245: [0.516246912003845,0.516246912217004]
                    246: [log ratio=,4.93705929488356]
                    247: [diag,18.6292,pfaffian,1.09207,ps,41.0026,rk,213.929]
                    248: 0.516246912217004
                    249: 266.4sec + gc : 8.277sec(276.8sec)
                    250: @end example
                    251:
                    252: @node n_wishartd.prob_by_ps,,, Numerical comptation of restricted function
                    253: @subsection @code{n_wishartd.prob_by_ps}
                    254: @findex n_wishartd.prob_by_ps
                    255:
                    256: @table @t
                    257: @item n_wishartd.prrob_by_ps(@var{m},@var{n},@var{[p1,p2,...]},@var{[s1,s2,...]},@var{t}[|@var{options}])
1.2       noro      258: computes the value of the distribution function of the largest eigenvalue of a Wishart matrix.
1.1       noro      259: @end table
                    260:
                    261: @table @var
1.2       noro      262: @item return
1.1       noro      263: @item m
1.2       noro      264: The number of variables.
1.1       noro      265: @item n
1.2       noro      266: The degrees of freedom.
                    267: @item [p1,p2,...,pk]
                    268: A list of the multiplicities of repeated eigenvalues.
                    269: @item [s1,s2,...,sk]
                    270: A list of repeated eigenvalues.
1.1       noro      271: @end table
                    272:
                    273: @itemize @bullet
                    274: @item
1.2       noro      275: This function compute a truncated power series solution up to a total degree
                    276: where the relative error between the current value and the previous value at the desired point
                    277: is less than @var{eps}. The default value of @var{eps} is @var{10^(-4)}.
                    278: The value of the distribution function is computed by using this power series.
                    279: @item If an option @var{eps} is set, the value is used as @var{eps}.
                    280: @item
                    281: If an option @var{eq} is not set,
                    282: a system of PDES satisfied by 1F1 on the diagonal region specified by @var{[p1,p2,...]} is computed.
                    283: If an option @var{eq} is set, the list specified by @var{eq} is regarded as the correct system of PDEs.
                    284: @item It is recommened to use this function when @var{t} is small.
1.1       noro      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,,, Numerical comptation of restricted function
                    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})
1.2       noro      305: computes a truncated power series solution up to the total degree @var{td}.
1.1       noro      306: @end table
                    307:
                    308: @table @var
                    309: @item return
1.2       noro      310: A list of polynomial
1.1       noro      311: @item z
1.2       noro      312: A list of differential operators with partial fraction coefficients.
1.1       noro      313: @item v
1.2       noro      314: A list of variables.
1.1       noro      315: @end table
                    316:
                    317: @itemize @bullet
                    318: @item
1.2       noro      319: The result is a list @var{[p,pd]} where
                    320: @var{p} is a truncated power series solution up to the total degree @var{td}
                    321: and @var{pd} is the @var{td} homogeneous part of @var{p}.
                    322: @item @var{z} cannot contain parameters other than the variables in @var{v}.
1.1       noro      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 Differential operators with partial fraction coefficients,,, n_wishartd.rr
                    339: @section Differential operators with partial fraction coefficients
                    340:
                    341: @menu
1.2       noro      342: * Representation of partial fractions::
                    343: * Representation of differential operators with partial fraction coefficients::
1.1       noro      344: * Operations on differential operators with partial fraction coefficients::
                    345: @end menu
                    346:
1.2       noro      347: @node Representation of partial fractions,,, Differential operators with partial fraction coefficients
                    348: @subsection Representation of partial fractions
1.1       noro      349:
1.2       noro      350: The coefficients of the PDE satisfied by the matrix 1F1 are written
                    351: as a sum of @var{1/yi} and @var{1/(yi-yj)} multiplied by constants.
                    352: Furthermore the result of diagonalization by l'Hopital rule
                    353: can also be written as a sum of partial fractions.
1.1       noro      354:
                    355: @itemize @bullet
                    356: @item
1.2       noro      357: A product @var{yi0^n0(yi1-yj1)^n1(yi2-yj2)^n2...(yik-yjk)^nk}
                    358: in the denominator of a fraction is represented as a list @var{[[yi0,n0],[yi1-yj1,n1],...,[yik-yjk,nk]]},
                    359: Where each @var{yi-yj} satisfies @var{i>j} and the factors are sorted according to an ordering.
1.1       noro      360: @item
1.2       noro      361: Let @var{f} be a power sum as above and @var{c} a constant.
                    362: Then a monomial @var{c/f} is represented by a list は @var{[c,f]}.
                    363: @var{f=[]} means that the denominator is 1.
1.1       noro      364: @item
1.2       noro      365: Finally @var{c1/f1+...+ck/fk} is represented as a list @var{[[c1,f1],...,[ck,fk]]},
                    366: where terms are sorted according to an ordering.
1.1       noro      367: @item
1.2       noro      368: We note that it is possible that a partial fraction is reduced to 0.
1.1       noro      369: @end itemize
                    370:
1.2       noro      371: @node Representation of differential operators with partial fraction coefficients,,, Differential operators with partial fraction coefficients
                    372: @subsection Representation of differential operators with partial fraction coefficients
1.1       noro      373:
1.2       noro      374: By using partial fractions explained in the previous section,
                    375: differential operators with partial fraction coefficients are represented.
                    376: Let @var{f1,...,fk} be partial fractions and @var{d1,...,dk} distributed monomials such that
                    377: @var{d1>...>dk}) with respected to the current monomial ordering.
                    378: Then a differential operator @var{f1*d1+...+fk*dk} is represented as a list @var{[f1,d1],...[fk,dk]]}.
1.1       noro      379:
                    380: @node Operations on differential operators with partial fraction coefficients,,, Differential operators with partial fraction coefficients
                    381: @subsection Operations on differential operators with partial fraction coefficients
                    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:
                    391: @node n_wishartd.wsetup,,, Operations on differential operators with partial fraction coefficients
                    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.2       noro      401: A natural number.
1.1       noro      402: @end table
                    403:
                    404: @itemize @bullet
1.2       noro      405: @item This function sets a @var{m}-variate computational enviroment. The variables are @var{y0,y1,...,ym} and @var{dy0,...,dym},
                    406: where @var{y0, dy0} are dummy variables for intermediate computation.
1.1       noro      407: @end itemize
                    408:
                    409: @node n_wishartd.addpf,,, Operations on differential operators with partial fraction coefficients
                    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.2       noro      418: A differential operator with partial fraction coefficients.
1.1       noro      419: @item p1, p2
1.2       noro      420: Differential operators with partial fraction coefficients.
1.1       noro      421: @end table
                    422:
                    423: @itemize @bullet
1.2       noro      424: @item This function computes the sum of differential operators @var{p1} and @var{p2}.
1.1       noro      425: @end itemize
                    426:
                    427: @node n_wishartd.mulcpf,,, Operations on differential operators with partial fraction coefficients
                    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.2       noro      436: A differential operator with partial fraction coefficients.
1.1       noro      437: @item c
1.2       noro      438: A partial fraction.
1.1       noro      439: @item p
1.2       noro      440: Differential operators with partial fraction coefficients.
1.1       noro      441: @end table
                    442:
                    443: @itemize @bullet
1.2       noro      444: @item This function computes the product of a partial fraction @var{c} and a differential operator @var{p}.
1.1       noro      445: @end itemize
                    446:
                    447: @node n_wishartd.mulpf,,, Operations on differential operators with partial fraction coefficients
                    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.2       noro      456: A differential operator with partial fraction coefficients.
1.1       noro      457: @item p1, p2
1.2       noro      458: Differential operators with partial fraction coefficients.
1.1       noro      459: @end table
                    460:
                    461: @itemize @bullet
1.2       noro      462: @item This function computes the product of differential operators @var{p1} and @var{p2}.
1.1       noro      463: @end itemize
                    464:
                    465: @node n_wishartd.muldpf,,, Operations on differential operators with partial fraction coefficients
                    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.2       noro      474: A differential operator with partial fraction coefficients.
1.1       noro      475: @item y
1.2       noro      476: A variable.
1.1       noro      477: @item p
1.2       noro      478: A differential operator with partial fraction coefficients.
1.1       noro      479: @end table
                    480:
                    481: @itemize @bullet
1.2       noro      482: @item
                    483: This function computes the product of the differential operator @var{dy} corresponding to a variable @var{y} and @var{p}.
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:
                    496: @node Experimental implementation of Runge-Kutta methods ,,, n_wishartd.rr
                    497: @section Experimental implementation of Runge-Kutta methods
                    498:
                    499: @menu
                    500: * rk_ratmat::
                    501: @end menu
                    502:
                    503: @node rk_ratmat,,, Experimental implementation of Runge-Kutta methods
                    504:
1.2       noro      505: In the function @code{n_wishartd.ps_by_hgm}, after computing the Pfaffian matrices for
                    506: the sytem of PDEs on a diagonal region, it executes a built-in function
                    507: @code{rk_ratmat} which computes an approximate solution of the Pfaffian system
                    508: by Runge-Kutta method for a spcified step size.
                    509: This function is repeated until the result gets stabilized, by doubling the step size.
                    510: @code{rk_ratmat} can be used as a general-purpose Runge-Kutta driver and we explain how to use it.
1.1       noro      511:
                    512: @subsection @code{rk_ratmat}
                    513: @findex rk_ratmat
                    514:
                    515: @table @t
                    516: @item rk_ratmat(@var{rk45},@var{num},@var{den},@var{x0},@var{x1},@var{s},@var{f0})
1.2       noro      517: solves a system of linear ODEs with rational function coefficients.
1.1       noro      518: @end table
                    519:
                    520: @table @var
                    521: @item return
1.2       noro      522: A list of real numbers.
1.1       noro      523: @item rk45
1.2       noro      524: 4 or 5.
1.1       noro      525: @item num
1.2       noro      526: An array of constant matrices.
1.1       noro      527: @item den
1.2       noro      528: A polynomial.
1.1       noro      529: @item x0, x1
1.2       noro      530: Real numbers.
1.1       noro      531: @item s
1.2       noro      532: A natural number.
1.1       noro      533: @item f0
1.2       noro      534: A real vector.
1.1       noro      535: @end table
                    536:
                    537: @itemize @bullet
                    538: @item
1.2       noro      539: Let @var{k} be the size of an array @var{num}.
                    540: The function @code{rk_ratmat} solves  an initial value problem
                    541: @var{dF/dx = P(x)F}, @var{F(x0)=f0} for @var{P(x)=1/den(num[0]+num[1]x+...+num[k-1]x^(k-1))} by a Runge-Kutta method.
1.1       noro      542: @item
1.2       noro      543: @var{rk45} specifies the order of a Runge-Kutta method. Adaptive methods are not implemented.
1.1       noro      544: @item
1.2       noro      545: The step size is specified by @var{s}. The step width is @var{(x1-x0)/s}.
1.1       noro      546: @item
1.2       noro      547: If the size of @var{f0} is @var{n},  each component of @var{num} is a square matrix of size @var{n}.
1.1       noro      548: @item
1.2       noro      549: The result is a list of real numbers @var{[r1,...,rs]} of length @var{s}.
                    550: @var{ri} is the 0-th component of the solution vector after the step @var{i}.
                    551: Before going to the next step the solution vector is divided by @var{ri}.
                    552: Therefore the 0-th component of the final solution vector [var{F(x1)} is equal to @var{rs*r(s-1)*...*r1}.
                    553: @item Since the ODE is linear, each step of Runge-Kutta method is also linear.
                    554: This enables us to apply a normalization such that the 0-th
                    555: component of each intermediate solution vector is set to 1.  By
                    556: applying this normalization we expect that all the components of
                    557: intermediate solution vectors can be represented by the format of
                    558: double precision floating point number.
1.3       noro      559: If there exist some components which cannot be represented by the format
                    560: of double precision floating point number in the initial vector @var{f0}, we apply this normalization
1.2       noro      561: to @var{f0}. After applying @code{rk_ratmat} we multiply the result for the normalized @var{f0} and the 0-th component
                    562: of the original @var{f0} to get the desired result.
1.1       noro      563: @end itemize
                    564:
                    565: @example
                    566: [...] F=ltov([sin(1/x),cos(1/x),sin(1/x^2),cos(1/x^2)]);
                    567: [ sin((1)/(x)) cos((1)/(x)) sin((1)/(x^2)) cos((1)/(x^2)) ]
                    568: [...] F0=map(eval,map(subst,F,x,1/10));
                    569: [ -0.54402111088937 -0.839071529076452 -0.506365641109759 0.862318872287684 ]
                    570: [...] N0=matrix(4,4,[[0,0,0,0],[0,0,0,0],[0,0,0,-2],[0,0,2,0]])$
                    571: [...] N1=matrix(4,4,[[0,-1,0,0],[1,0,0,0],[0,0,0,0],[0,0,0,0]])$
                    572: [...] N=ltov([N0,N1])$
                    573: [...] D=x^3$
                    574: [...] R=rk_ratmat(5,N,D,1/10,10,10^4,F0)$
                    575: [...] for(T=R,A=1;T!=[];T=cdr(T))A *=car(T)[1];
                    576: [...] A;
                    577: 0.0998334
                    578: [...] F1=map(eval,map(subst,F,x,10));
                    579: [ 0.0998334166468282 0.995004165278026 0.00999983333416666 0.999950000416665 ]
                    580: @end example
                    581:
                    582:
                    583: @comment --- おまじない ---
                    584: @node Index,,, Top
                    585: @unnumbered Index
                    586: @printindex fn
                    587: @printindex cp
                    588: @iftex
                    589: @vfill @eject
                    590: @end iftex
                    591: @summarycontents
                    592: @contents
                    593: @bye
                    594: @comment --- おまじない終り ---
                    595:

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