home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / INTEGRAT.PAK / LOGCASES.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  6.2 KB  |  194 lines
  1.  
  2. (****************************************************************************
  3. *               Evaluation Multiple Poles Meijer's G-function
  4. *
  5. *****************************************************************************)
  6.  
  7. (* HISTORY:
  8.     This code used for evaluation poles Meijer's G-function in the 
  9.     following cases:
  10.     a) finite number of multiple poles any order;
  11.     b) infinite number of first order poles.
  12.  
  13.     In October 1991 was added some code to evaluate 
  14.     c) infinite number of second order poles.
  15.  
  16. *)
  17.  
  18.  
  19.  MeijerLogCase[parn_,parp_,parm_,parq_,z_ ] :=
  20.   Module[ {answer={},grpoles,grzero,i}, 
  21.      grpoles = FindGroupsPoles[ parm ];
  22.      For[ i=1,i<=Length[grpoles],i++,
  23.           grzero = OneGroupZero[ grpoles[[i,1]],parp,{} ][[3]];
  24.           If[ Length[grzero] != 0,
  25.               answer = Append[answer,OrderPoles[
  26.                               Join[ Map[ {#,1}&, grpoles[[i]] ],
  27.                                     Map[ {#,-1}&,grzero       ]]]],
  28.               answer = Append[answer,OrderPoles[
  29.                                     Map[ {#,1}&, grpoles[[i]] ]]]
  30.             ]
  31.         ];
  32.   answer = 
  33.   Apply[ Plus,Map[ 
  34.    EvalResidues[#,parn,parp,parm,parq,z/.arg->Arg]&,CountMultiple[answer] ] ];
  35.   answer
  36.     ]       
  37.  
  38.  
  39.  FindGroupsPoles[v_] := OneGroup[ v[[1]], Rest[v], {v[[1]]} ]
  40.  
  41.  FindGroupsPoles[v_] := {} /; Length[v] == 0
  42.  
  43.  OneGroup[ a_,{v1___,b_,v2___},c_ ] :=
  44.      OneGroup[ a, {v1,v2}, Append[c,b] ] /; IntegerQ[Expand[a-b]]
  45.  
  46.  OneGroup[ a_,b_,c_ ] := Append[ FindGroupsPoles[b],c ] 
  47.  
  48.  OneGroupZero[ a_,{v1___,b_,v2___},c_ ] :=
  49.      OneGroupZero[ a, {v1,v2}, Append[c,b] ] /; IntegerQ[Expand[a-b]]
  50.    
  51.    
  52.  OrderPoles[{c1___,a_,c2___,b_,c3___}] :=
  53.                OrderPoles[{c1,b,c2,a,c3}] /; Expand[a[[1]]-b[[1]]] > 0
  54.  
  55.  OrderPoles[{c___}] := {c}
  56.  
  57.  CountMultiple[ v_ ] := Flatten[ Map[CountInGroup[#]&,v],1 ]
  58.  
  59.  
  60.  CountInGroup[v_] :=
  61.   Module[ {u=v,kr=0,answer={}},
  62.        While[ Length[u] !=0,
  63.           If[ u[[1,2]] == 1, kr += 1, kr -= 1 ];
  64.           If[ kr > 0,
  65.              If[Length[Rest[u]] != 0,
  66.                 PrependTo[answer,{kr,Expand[u[[2,1]]-u[[1,1]]],-u[[1,1]]}],
  67.         If[ kr == 1, 
  68.             PrependTo[ answer, {kr,"Infinity",u[[1,1]]} ],
  69.            If[ kr == 2,
  70.                PrependTo[answer,{kr,"Infinity",u[[1,1]],
  71.                  answer[[1,2]]} ],
  72.                 answer = FailInt; u = {FailInt}
  73.            ]
  74.         ]
  75.          ]
  76.      ];
  77.        u = Rest[u] ];
  78.        answer
  79.  ]
  80.  
  81.  
  82.  EvalResidues[{1,"Infinity",at_},parn_,parp_,parm_,parq_,z_] := 
  83.    If[MemberQ[parm,at],
  84.       FinalGfunToHyper1[
  85.            at,parn,parp,First[ReducePar[parm,{at}]],parq,z],
  86.       SpecialTransf[at,
  87.            parn,Append[parp,at],Append[parm,at],parq,z,1]
  88.      ]  
  89.  
  90.  FinalGfunToHyper1[ at_, n_,{v3___,a_,v4___},{v1___,b_,v2___},q_,z_ ] :=
  91.      (-1)^(a-b) *
  92.      FinalGfunToHyper[ at, Join[n, {a}],{v3,v4},{v1,v2},Join[q,{b}],z ] /;
  93.  IntegerQ[a-b] && Positive[a-b]
  94.  
  95.  FinalGfunToHyper1[ w__ ] := FinalGfunToHyper[w]
  96.   
  97.  EvalResidues[{2,"Infinity",at_,m_},parn_,parp_,parm_,parq_,z_] := 
  98.    If[ MemberQ[parm,at],
  99.        SecondOrder[ at,m,Complement[parm,{at,at-m}],1-parn,parp,1-parq,z],
  100.        SpecialTransf[at,parn,Join[parp,{at,at}],
  101.             Join[parm,{at,at}],parq,z,2]
  102.    ]
  103.  
  104.  SpecialTransf[at_,parn_,{v1___,b_,v2___},{v3___,a_,v4___},parq_,z_,k_] :=
  105.    (-1)^(b-a) SpecialTransf[at, Append[parn,b],{v1,v2},{v3,v4},
  106.                              Append[parq,a],z,k ] /;
  107.   IntegerQ[Expand[b-a]] && Expand[b-a] > 0
  108.  
  109.      
  110.  SpecialTransf[at_,parn_,parp_,parm_,parq_,z_,k_] :=
  111.    Module[ {r},
  112.      r = ReducePar[parn,parq];
  113.      If[ k==1,
  114.      FinalGfunToHyper[
  115.         at,r[[1]],parp,First[ReducePar[parm,{at}]],r[[2]],z],
  116.      SecondOrder[ at,0,Complement[parm,{at,at}],
  117.             1-r[[1]],parp,1-r[[2]],z],
  118.        ]
  119.    ]
  120.  
  121.  EvalResidues[{order_,0,at_},__] := 0
  122.  
  123.  EvalResidues[{order_,1,at_},parn_,parp_,parm_,parq_,z_] := 
  124.   Module[ {rr,eps},
  125.    rr =  
  126.      Apply[ Times,
  127.            Map[GammaRes[#,-eps,order]&,1-parn-at] ] *
  128.      Apply[ Times,
  129.            Map[GammaRes[#, eps,order]&, parm+at ] ] *
  130.      Apply[ Times,
  131.            Map[GammaResInv[#, eps,order]&, parp+at ] ] *
  132.      Apply[ Times,
  133.            Map[GammaResInv[#,-eps,order]&,1-parq-at] ] *
  134.      Series[Exp[-eps Log[z]],{eps,0,order-1}] z^(-at);
  135.   Coefficient[ rr,eps,Exponent[rr,eps] ]
  136.   ] 
  137.  
  138.  EvalResidues[{order_,k_Integer/;k>1,at_},parn_,parp_,parm_,parq_,z_] :=
  139.     Sum[ EvalResidues[{order,1,at-j},parn,parp,parm,parq,z],
  140.      {j,0,k-1} ]
  141.  
  142.  GammaRes[u_Integer/;u<=0,e_,1] := 
  143.    SeriesData[e, 0, {1/((-1)^u*Gamma[1 - u])}, 0, 1, 1]  
  144.  
  145.  GammaRes[u_Integer/;u<=0,e_,2] := SeriesData[e, 0, {1/((-1)^u*Gamma[1 - u]), 
  146.     PolyGamma[0, 1 - u]/((-1)^u*Gamma[1 - u])}, 0, 2, 1]
  147.  
  148.  GammaRes[0,e_,n_] := GammaRes[1,e,n]
  149.  
  150.  GammaRes[u_Integer/;u<0,e_,n_] := 
  151.     (-1)^(-u) GammaRes[1,e,n] GammaRes[1,-e,n] GammaResInv[1-u,-e,n] 
  152.  
  153.  GammaRes[u_,e_,n_] := Series[Gamma[u+e],
  154.                              {If[Head[e] === Symbol,e,-e],0,n-1}]  
  155.  
  156.  GammaResInv[u_,e_,n_] := Series[1/Gamma[u+e],
  157.                              {If[Head[e] === Symbol,e,-e],0,n-1}]  
  158.  
  159.  SecondOrder[ at_, m_, a_, b_, c_, d_, z_ ] :=
  160.   Module[ { const },
  161.     const = z^at (-1)^m *
  162.         SimplifyGamma[Times@@(MultGamma[#]&/@{a-at,b+at})/(
  163.                   Times@@(MultGamma[#]&/@{c-at,d+at}) m!)];
  164.     -const Log[z] *
  165.     HyperInteg[ Join[b+at,1+at-c],Join[d+at,1+at-a,{m+1}],
  166.                (-1)^(Length[c]-Length[a]) z ] 
  167.     +
  168.     const *
  169.     sum@@{ ((-1)^If[EvenQ[Length[c]-Length[a]],0,1] z)^dummy/dummy! *
  170.      Times@@(MultPochham[#,dummy]&/@{b+at,1+at-c})/(
  171.      Times@@(MultPochham[#,dummy]&/@{d+at,1+at-a,{m+1}}) ) *
  172.          Expand[ SimplifyPolyGammaS[dummy,
  173.       Plus@@(PolyGamma[0,#]&/@Expand[d+at+dummy]) - 
  174.       Plus@@(PolyGamma[0,#]&/@Expand[b+at+dummy]) - 
  175.       Plus@@(PolyGamma[0,#]&/@Expand[c-at-dummy]) + 
  176.       Plus@@(PolyGamma[0,#]&/@Expand[a-at-dummy]) + 
  177.       PolyGamma[0,1+m+dummy] + PolyGamma[0,1+dummy] ] ],
  178.     {dummy, 0, Infinity}}
  179.  ]
  180.  
  181.  SimplifyPolyGammaS[ l_, expr_ ] :=
  182.    Module[ { r = SimplifyPolyGamma[expr] },
  183.       If[ FreeQ[r,PolyGamma[0,a_. + b_?Negative l]],
  184.       r,
  185.       r/.PolyGamma[0,a_. + b_?Negative l]:>
  186.       PolyGamma[0,-a-b l]-1/Factor[a+b l]-Pi Cot[Expand[Pi (a+b l)]]
  187.       ]/.
  188.     { Cot[r_ + p_ l] :> Cot[r]/; FreeQ[r,l] && Abs[p]==Pi,
  189.       Cot[p_ l] :> 0,
  190.       Tan[r_. + p_ l] :> Tan[r]/; FreeQ[r,l] && Abs[p]==Pi
  191.         }
  192.    ]
  193.  
  194.