# to unbundle, sh this file (in an empty directory)
# or use ftp://netlib.bell-labs.com/netlib/access/unshar.c.gz
echo mathematica 1>&2
sed >mathematica <<'-------cut here----- mathematica' 's/^X//'
X(*:Name: `RknStabInt` *)
X
X(*:Title:    Computation of the interval of stability of
X             Runge-Kutta-Nystrom methods *)
X
X(*:Authors:  Beatrice Paternoster,
X             Dipartimento di Informatica e Applicazioni,
X             Universita' di Salerno, Italy
X             E-mail: beapat@dia.unisa.it
X
X             Massimo Cafaro,
X             Facolta' di Ingegneria,
X             Universita' di Lecce, Italy
X             E-mail: cafaro@sara.unile.it *)
X
X(*:Keywords: A-stability, P-stability, zero dissipative methods *)
X
X(*:Requirements: none. *)
X
X(*:Warnings: filenames can't contain spaces*)
X
X(*:Sources:   P.J. Van Der Houwen,B.P. Sommeijer,Nguyen Huu Cong,
X              Stability of collocation-based RKN methods,
X              BIT 31  (1991) 469--481
X
X              B.Paternoster,M.Cafaro,
X              Computation of the interval of stability of
X              Runge-Kutta-Nystrom methods, to appear in J.Symb.Comp.*)
X
XBeginPackage["`RknStabInt`"]
X
XStabInterval::usage = "StabInterval[filename] return the stability
X             interval of the RKN method saved in the file filename;
X             the RKN coefficients must be saved as matrices in the
X             following order: A b bp c"
X
XBegin["`Private`"]
X
X
X(* ReadData reads the RKN coefficients from the file filename *)
X
XReadData[filename_]:=Module[{name},
X         name=OpenRead[ToString[filename]];
X         A=Read[name,Expression];
X         b=Read[name,Expression];
X         bp=Read[name,Expression];
X         c=Read[name,Expression];
X         Close[name]];
X
X
X
X(* BuildMatrix computes the stability matrix M of the RKN method,
X   the determinant det and the trace *)
X
XBuildMatrix[A_,b_,bp_,c_]:=Module[
X         {stages, e, m11, m12, m21, m22, M},
X         stages=Length[A];
X         inv=Chop[Inverse[IdentityMatrix[stages]-z A]];
X         e=Transpose[{Table[i-i+1,{i,stages}]}];
X         m11=1+ z ((b . inv . e) [[1,1]]);
X         m12=1+ z ((b . inv . Transpose[c]) [[1,1]]);
X         m21=z ((bp . inv . e) [[1,1]]);
X         m22=1+ z ((bp . inv . Transpose[c]) [[1,1]]);
X         M={{ m11,m12 },{ m21,m22 }};
X         det=Together[Det[M]];
X         trace=Together[Sum[M[[i,i]],{i,2}]]];
X
X
X
X(* BuildExpr determines if the RKN method is zero-dissipative,
X   and constructs the three expressions necessary to compute
X   the stability interval *)
X
XBuildExpr[d_,t_]:=Module[{},
X         If[!IntegerQ[d],
X         numexp=3;
X         exp[1]=Together[(t)^2-((d)+1)^2];
X         exp[2]=Together[(d)-1];
X         exp[3]=Together[(-d)-1],
X         numexp=1;
X         Print["zero dissipative method"];
X         exp[1]=Together[(t)^2-(4(d))]]];
X
X
X
X(* MakeRoots returns a list containing the roots of det and trace *)
X
XMakeRoots[exp_]:=Module[{nu,de,root,zero},
X         nu=Numerator[exp];
X         de=Denominator[exp];
X         zero=Flatten[z /. NSolve[nu==0,z]];
X         If[!IntegerQ[de],
X         root=Flatten[z /. NSolve[de==0,z]],
X         root={}];
X         temp=N[Union[root,zero]]];
X
X
X(* NegativeRoots extracts the negative roots from the list, and puts
X   them in a list of decreasing order starting from zero *)
X
XNegativeRoots[list_]:=Module[{neg,numzer},
X         neg={0};
X         roots=Reverse[N[Map[Chop,list]]];
X         numzer=Length[roots];
X         Do[If[Sign[roots[[i]]]==-1,
X         AppendTo[neg,roots[[i]]]],{i,numzer}];
X         roots=neg];
X
X
X
X(* Bound computes the largest interval of the negative real axis, in which
X   |det| <= 1 and trace^2-4 det < 0, using the list of negative roots
X   constructed by NegativeRoots. If this list is empty, Bound decides
X   if the interval of stability is empty or non limited *)
X
XBound[exp_]:=Module[{numzer},
X         Clear[z];
X         numzer=Length[roots];
X         If[numzer==1,EmptyList[exp],NotEmptyList[numzer,exp]];
X         Clear[z]];
X
X         EmptyList[exp_]:=Module[{},z=-10;
X         If[Sign[N[exp]]==-1,res=-Infinity,res=0]];
X
X         NotEmptyList[nz_,exp_]:=Module[{i,negative},
X         i=2;negative=True;
X         While[i<=nz && negative==True,
X         ComputeSign[i,exp];
X         If[s==-1,Clear[z];i++,negative=False]];
X         If[i==nz+1,Clear[z];CheckLastRoot[nz,exp],
X         res=roots[[i-1]]]];
X
X         ComputeSign[k_,exp_]:=Module[{},
X         z=(roots[[k]]+roots[[k-1]])/2;
X         s=Sign[N[exp]]];
X
X         CheckLastRoot[nz_,exp_]:=Module[{},
X         z=roots[[nz]]-10;
X         If[Sign[N[exp]]==-1,res=-Infinity,res=roots[[nz]]]];
X
X(* StabInterval determines the interval of stability of the RKN method,
X   intersecting the interval computed by Bound *)
X
XStabInterval[filename_]:=Module[{stabint},
X         ReadData[filename];
X         BuildMatrix[A,b,bp,c];
X         BuildExpr[det,trace];
X         Do[MakeRoots[exp[i]];
X         NegativeRoots[temp];
X         Bound[exp[i]];
X         ext[i]=res,{i,numexp}];
X         stabint={Max[Flatten[Table[ext[i],{i,numexp}]]],0};
X         Print[StringForm["stability interval = ``",Expand[stabint]]];
X         ClearAll[A,b,bp,c,det,trace];
X         ClearAll[numexp,temp,roots,s,res]];
X
XEnd[]
X
XProtect[StabInterval]
X
XEndPackage[]
X
X
-------cut here----- mathematica
echo data1 1>&2
sed >data1 <<'-------cut here----- data1' 's/^X//'
X{{1/12-Sqrt[15]/60,0},{Sqrt[15]/60,1/12-Sqrt[15]/60}}
X{{0,1/2}}
X{{0,1}}
X{{1/2,1/2}}
-------cut here----- data1
echo data2 1>&2
sed >data2 <<'-------cut here----- data2' 's/^X//'
X{{0.052320267566927,0,0},
X{-0.17329232352333,0.052320267566927,0},
X{-0.01271397498318,0.043727040749588,0.052320267566927}}
X{{0,0,1/2}}
X{{0,0,1}}
X{{1/2,3/10,1/2}}
-------cut here----- data2
echo data3 1>&2
sed >data3 <<'-------cut here----- data3' 's/^X//'
X{{0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
X   {2.0 10^(-4),0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
X   {2.66666666666666666666666666667 10^(-4),
X    5.33333333333333333333333333333 10^(-4),
X    0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
X   {2.91666666666666666666666666667 10^(-3),
X    -4.16666666666666666666666666667 10^(-3),
X    6.25 10^(-3),0,0,0,0,0,0,0,0,0,0,0,0,0,0},
X   {1.64609053497942386831275720165 10^(-3),0,
X    5.48696844993141289437585733882 10^(-3),
X    1.75582990397805212620027434842 10^(-3),
X    0,0,0,0,0,0,0,0,0,0,0,0,0},
X   {1.9456 10^(-3),0,
X    7.15174603174603174603174603175 10^(-3),
X    2.91271111111111111111111111111 10^(-3),
X    7.89942857142857142857142857143 10^(-4),
X    0,0,0,0,0,0,0,0,0,0,0,0},
X   {5.6640625 10^(-4),0,
X    8.80973048941798941798941798942 10^(-4),
X    -4.36921296296296296296296296296 10^(-4),
X    3.39006696428571428571428571429 10^(-4),
X    -9.94646990740740740740740740741 10^(-5),
X    0,0,0,0,0,0,0,0,0,0,0},
X   {3.08333333333333333333333333333 10^(-3),0,0,
X    1.77777777777777777777777777778 10^(-3),2.7 10^(-3),
X    1.57828282828282828282828282828 10^(-3),
X    1.08606060606060606060606060606 10^(-2),
X    0,0,0,0,0,0,0,0,0,0},
X   {3.65183937480112971375119150338 10^(-3),0,
X    3.96517171407234306617557289807 10^(-3),
X    3.19725826293062822350093426091 10^(-3),
X    8.22146730685543536968701883401 10^(-3),
X    -1.31309269595723798362013884863 10^(-3),
X    9.77158696806486781562609494147 10^(-3),
X    3.75576906923283379487932641079 10^(-3),
X    0,0,0,0,0,0,0,0,0},
X   {3.70724106871850081019565530521 10^(-3),0,
X    5.08204585455528598076108163479 10^(-3),
X    1.17470800217541204473569104943 10^(-3),
X    -2.11476299151269914996229766362 10^(-2),
X    6.01046369810788081222573525136 10^(-2),
X    2.01057347685061881846748708777 10^(-2),
X    -2.83507501229335808430366774368 10^(-2),
X    1.48795689185819327555905582479 10^(-2),
X    0,0,0,0,0,0,0,0},
X   {3.51253765607334415311308293052 10^(-2),0,
X    -8.61574919513847910340576078545 10^(-3),
X    -5.79144805100791652167632252471 10^(-3),
X    1.94555482378261584239438810411,
X    -3.43512386745651359636787167574,
X    -0.109307011074752217583892572001,
X    2.3496383118995166394320161088,
X    -0.756009408687022978027190729778,
X    0.109528972221569264246502018618,0,0,0,0,0,0,0},
X   {2.05277925374824966509720571672 10^(-2),0,
X    -7.28644676448017991778247943149 10^(-3),
X    -2.11535560796184024069259562549 10^(-3),
X    0.927580796872352224256768033235,
X    -1.65228248442573667907302673325,
X    -2.10795630056865698191914366913 10^(-2),
X    1.20653643262078715447708832536,
X    -0.413714477001066141324662463645,
X    9.07987398280965375956795739516 10^(-2),
X    5.35555260053398504916870658215 10^(-3),0,0,0,0,0,0},
X   {-0.143240788755455150458921091632,0,
X    1.25287037730918172778464480231 10^(-2),
X    6.82601916396982712868112411737 10^(-3),
X    -4.79955539557438726550216254291,
X    5.69862504395194143379169794156,
X    0.755343036952364522249444028716,
X    -0.127554878582810837175400796542,
X    -1.96059260511173843289133255423,
X    0.918560905663526240976234285341,
X    -0.238800855052844310534827013402,
X    0.159110813572342155138740170963,0,0,0,0,0},
X   {0.804501920552048948697230778134,0,
X    -1.66585270670112451778516268261 10^(-2),
X    -2.1415834042629734811731437191 10^(-2),
X     16.8272359289624658702009353564,
X    -11.1728353571760979267882984241,
X    -3.37715929722632374148856475521,
X    -15.2433266553608456461817682939,
X    17.1798357382154165620247684026,
X    -5.43771923982399464535413738556,
X    1.38786716183646557551256778839,
X    -0.592582773265281165347677029181,
X    2.96038731712973527961592794552 10^(-2),0,0,0,0},
X   {-0.913296766697358082096250482648,0,
X    2.41127257578051783924489946102 10^(-3),
X    1.76581226938617419820698839226 10^(-2),
X    -14.8516497797203838246128557088,
X    2.15897086700457560030782161561,
X    3.99791558311787990115282754337,
X    28.4341518002322318984542514988,
X    -25.2593643549415984378843352235,
X    7.7338785423622373655340014114,
X    -1.8913028948478674610382580129,
X    1.00148450702247178036685959248,
X    4.64119959910905190510518247052 10^(-3),
X    1.12187550221489570339750499063 10^(-2),0,0,0},
X   {-0.275196297205593938206065227039,0,
X    3.66118887791549201342293285553 10^(-2),
X    9.7895196882315626246509967162 10^(-3),
X    -12.293062345886210304214726509,
X    14.2072264539379026942929665966,
X    1.58664769067895368322481964272,
X    2.45777353275959454390324346975,
X    -8.93519369440327190552259086374,
X    4.37367273161340694839327077512,
X    -1.83471817654494916304344410264,
X    1.15920852890614912078083198373,
X    -1.72902531653839221518003422953 10^(-2),
X    1.93259779044607666727649875324 10^(-2),
X    5.20444293755499311184926401526 10^(-3),0,0},
X   {1.30763918474040575879994562983,0,
X    1.73641091897458418670879991296 10^(-2),
X    -1.8544456454265795024362115588 10^(-2),
X    14.8115220328677268968478356223,
X    9.38317630848247090787922177126,
X    -5.2284261999445422541474024553,
X    -48.9512805258476508040093482743,
X    38.2970960343379225625583875836,
X    -10.5873813369759797091619037505,
X    2.43323043762262763585119618787,
X    -1.04534060425754442848652456513,
X    7.17732095086725945198184857508 10^(-2),
X    2.16221097080827826905505320027 10^(-3),
X    7.00959575960251423699282781988 10^(-3),0,0}}
X{{1.21278685171854149768890395495 10^(-2),0,0,0,0,0,
X    8.62974625156887444363792274411 10^(-2),
X    0.252546958118714719432343449316,
X    -0.197418679932682303358307954886,
X    0.203186919078972590809261561009,
X    -0.0207758080777149166121933554691,
X    0.109678048745020136250111237823,
X    0.0380651325264665057344878719105,
X    0.0116340688043242296440927709215,
X    4.65802970402487868693615238455 10^(-3),0,0}}
X{{1.21278685171854149768890395495 10^(-2),0,0,0,0,0,
X     0.0908394342270407836172412920433,
X     0.315683697648393399290429311645,
X     -0.263224906576909737811077273181,
X     0.304780378618458886213892341513,
X     -0.415516161554298332243867109382,
X     0.246775609676295306562750285101,
X     0.152260530105866022937951487642,
X     0.0814384816302696075086493964505,
X     0.0850257119389081128008018326881,
X     -9.15518963007796287314100251351 10^(-3),0.025}}
X{{0,0.02,0.04,0.1,
X    0.133333333333333333333333333333,0.16,0.05,0.2,0.25,
X    0.333333333333333333333333333333,0.5,
X    0.555555555555555555555555555556,0.75,
X    0.857142857142857142857142857143,
X    0.945216222272014340129957427739,1.0,1.0}}
-------cut here----- data3