(* Table of singularities with their weights for quadratic and differential approximants *)
(* Options *)
(*! was 15 *)mindlist = 20; (* number of approximants *)
ncrop = 22; (* disregard coefficients with N>ncrop *)
ncrop = 38; (* disregard coefficients with N>ncrop *)
ndrop = 0; (* how many coeff. to drop from beginning of series *)
xyrange = 2.5/asave[name]; (* plot range *)
sect = "sngtable";
entity = "table";
cf = If[nm > ncrop, Take[func, ncrop], func];
cf = Drop[cf, ndrop];
nmax = Length[cf];
(* Create table of indexes of approximants *)
indlist = {};
Do[Do[
ind = {n1, n2, n3};
nuse = Plus@@ind + 2;
prior = 1.1^nuse;
If[n1 < 0 || n2 < 0 || n3 < 0 || nuse > nmax,
Continue[]];
If[Abs[n1 - n2] > 1 || Abs[n2 - n3] > 1 || Abs[n3 - n1] > 1,
Continue[]];
If[!(n3<=n2<=n1<=n3+1), Continue[]];
indlist = Append[indlist, {ind, prior, nuse}];
, {n2, 0, n1}, {n3, 0, n2}], {n1, 0, nm}];
indlist = Sort[indlist, #1[[2]] > #2[[2]] &];
mi = Length[indlist];
If[mindlist < mi, indlist = Take[indlist, mindlist]; mind = mindlist, mind = mi];
indlist = Sort[indlist, #1[[3]] < #2[[3]] &];
singsf = ToFileName[dir, "sings.htm"];
of = open[singsf];
"<!DOCTYPE HTML PUBLIC \"-//W3C//DTD HTML 3.2//EN\">\n
<HTML><HEAD>\n
<TITLE>Singularities of Moller-Plesset series: example "" <> name <>
""</TITLE>\n
<META name=\"description\" content=\"Singularities of approximants for Moller-Plesset perturbation series\">\n
<META name=\"keywords\" content=\"quantum mechanics, Moller Plesset perturbation theory, summation, Pade approximants, branch cuts\">\n
</HEAD><BODY BGCOLOR=\"#C0C0C0\"><A NAME=\"top\"></A>\n
<H1>Singularities of Møller-Plesset series: example "" <> name <>
""</H1>\n
<H2>Molecule " <> mol <> ". Basis "<> bas <> ". Structure "" <> fdat <>
""</H2>" // p;
"<H3>Content</H3>\n
<UL>\n
<LI><A HREF=\"#quadr\">Definition of quadratic approximants</A>\n
<UL>" // p;
Clear[pol, c, funs];
Do[{{n1, n2, n3}, prior, nc} = indlist[[nind]];
indtxt = "[" <> ToString[n1] <> ", " <> ToString[n2] <> ", " <> ToString[n3] <> "]";
"<LI><A HREF=\"#quadr" <> ToString[nind] <> "\">Approximant " <> indtxt <> "</A>" // p;
c[0] = 1;
pol[0] = Sum[c[n]z^n, {n, 0, n1}];
pol[1] = Sum[c[n + n1 + 1]z^n, {n, 0, n2}];
pol[2] = Sum[c[n + n1 + n2 + 2]z^n, {n, 0, n3}];
funs[z_] = Sum[cf[[n]]z^(n - 1), {n, nc}];
(*Quadratic approximant*)
pattern = pol[0] + pol[1]f[z] + pol[2]f[z]^2;
eqs = Take[CoefficientList[pattern /. f -> funs, z], nc];
s = Solve[eqs == Table[0, {nc}], Table[c[n], {n, nc}]];
ls = Length[s];
If[ls == 0, Print["No solutions of linear equations found!"]];
If[ls > 1, Print["More than one solution of linear equations!"]];
s1 = s[[1]];
{pa, pb, pc} = {pol[2], pol[1], pol[0]} /. s1;
pold = pb^2 - 4 pa pc;
sn = NSolve[pold == 0, z];
sn = Sort[sn, Abs[((z /. #1) - 0)] < Abs[((z /. #2) - 0)] &];
sqposit[nind] = sings = z /. sn;
If[nc<=5, singularities[nser,nc]=sings//N; (*Print["Sing-s(",nc,"): ",singularities[nser,nc]]*)];
sqweight[nind] = 1/2 Sqrt[-z D[pold/pa^2, z]] /. sn;
, {nind, mind}];
"</UL>" // p;
"</UL>" // p;
rulersing = StringReplace[ruler, "/index.htm\">"->"/sings.htm\">"];
"<HR>" <> rulersing <> "<HR>" // p;
"<TABLE><TR ALIGN=\"CENTER\">\n
<TD ALIGN=\"CENTER\" VALIGN=\"CENTER\"><A HREF=\"index.htm#sings3\">Plot of singularities</A></TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"TOP\"><IMG SRC=\"../../../../blank.gif\" WIDTH=32 HEIGHT=43 ALT=\"Blank\"></TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"CENTER\"><A HREF=\"../index.htm\">" <>
"<IMG SRC=\"../../../../iconz/molec1.gif\" WIDTH=36 HEIGHT=41 BORDER=0 ALIGN=\"ABSMIDDLE\" ALT=\"Molecule - icon for Allen-data\">" <>
"List of examples</A></TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"TOP\"><IMG SRC=\"../../../../blank.gif\" WIDTH=32 HEIGHT=43 ALT=\"Blank\"></TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"CENTER\"><A HREF=\"../proglist.htm\"><I>Mathematica</I> programs</A></TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"TOP\"><IMG SRC=\"../../../../blank.gif\" WIDTH=32 HEIGHT=43 ALT=\"Blank\"></TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"CENTER\"><A HREF=\"../../index.htm\">Work in UMassD</A></TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"TOP\"><IMG SRC=\"../../../../blank.gif\" WIDTH=32 HEIGHT=43 ALT=\"Blank\">
</TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"CENTER\"><A HREF=\"../../../reports.htm\">" <>
"<IMG SRC=\"../../../waste.gif\" WIDTH=32 HEIGHT=43 BORDER=0 ALIGN=\"ABSMIDDLE\" ALT=\"Waste icon\">Unpublished reports</A></TD>\n
</TR></TABLE><HR>" // p;
(* Definition of quadratic approximant *)
"<A NAME=\"quadr\"></A><H3>Quadratic approximants</H3>\n
<P>[<I>n</I><SUB>1</SUB>, <I>n</I><SUB>2</SUB>, <I>n</I><SUB>3</SUB>] approximant is defined\n
as a solution of the quadratic equation\n
<BR><I>A</I><SUB></SUB>(<I>z</I>)<I>f</I><SUP>2</SUP> +
<I>B</I>(<I>z</I>)<I>f</I> +
<I>C</I>(<I>z</I>) = 0\n
<BR>with polynomial coefficients <I>A</I>(<I>z</I>),\n
<I>B</I>(<I>z</I>) and\n
<I>C</I>(<I>z</I>) of degree\n
<I>n</I><SUB>3</SUB>, <I>n</I><SUB>2</SUB> and <I>n</I><SUB>1</SUB> respectively.\n
<P>Square-root singularities are determined as zeroes of the discriminant\n
<BR><I>D</I>(<I>z</I>) = <I>B</I><SUP>2</SUP>(<I>z</I>) - 4<I>A</I>(<I>z</I>)<I>C</I>(<I>z</I>).\n
<BR>The weight <I>c</I> of the singularity <I>z</I><SUB>c</SUB> is defined so that\n
<BR><I>f</I> ~ <I>c</I>(1 - <I>z</I>/<I>z</I><SUB>c</SUB>)<SUP>1/2</SUP>
at <I>z</I> -> <I>z</I><SUB>c</SUB>.\n
<BR>The weight is calculated by formula\n
<BR><I>c</I> = 1/2[-<I>z</I>(<I>D</I>/<I>A</I><SUP>2</SUP>)']<SUP>1/2</SUP>\n
<BR>where r. h. s. of the above equation is evaluated at <I>z</I> = <I>z</I><SUB>c</SUB>.
<BR>" // p;
Clear[pol, c, funs];
Do[{{n1, n2, n3}, prior, nc} = indlist[[nind]];
indtxt = "[" <> ToString[n1] <> ", " <> ToString[n2] <> ", " <> ToString[n3] <> "]";
"<BR><A NAME=\"quadr" <> ToString[nind] <> "\"></A><TABLE BORDER=0 CELLSPACING=1 CELLPADDING=5 WIDTH=\"780\">\n
<TR><TH ALIGN=\"CENTER\" VALIGN=\"MIDDLE\" BGCOLOR=\"#FFFFFF\" COLSPAN=4>
Table " <> ToString[nind] <> ". Singularities with their weights for the quadratic approximant "
<> indtxt <> "<BR>The most <FONT COLOR=\"#FF0000\">stable</FONT> singularity is highlighted.</TH></TR>\n
<TR ALIGN=\"CENTER\" VALIGN=\"MIDDLE\" BGCOLOR=\"#FEF8CB\"><TH WIDTH=\"40\">No.</TH>
<TH WIDTH=\"200\"><I>z</I><SUB>c</SUB></TH><TH WIDTH=\"200\"><I>c</I></TH>
<TH WIDTH=\"340\">Location in the complex plane</TH></TR>" // p;
sings = sqposit[nind]//Chop//N;
msing = Length[sings];
weights = sqweight[nind]//Chop//N;
(*Stable singularity*)
singsprev = If[nind == 1, sings, sqposit[nind - 1] // Chop // N];
singsnext = If[nind == mind, sings, sqposit[nind + 1] // Chop // N];
{mprev, mnext} = Length /@ {singsprev, singsnext};
stable = {Null, Null, Infinity};
Do[sing = sings[[n]];
If[Im[sing] < 0, Continue[]];
err1 = Abs[singsprev[[nprev]] - sing] + Abs[singsnext[[nnext]] - sing];
If[err1 < stable[[3]], stable = {sing, n, err1}], {n, msing}, {nprev,
mprev}, {nnext, mnext}];
{sing0, nsing0, err0} = stable;
Do[
hstart = If[nsing==nsing0, "<B>", ""];
hend = If[nsing==nsing0, "</B>", ""];
"<TR ALIGN=\"CENTER\" VALIGN=\"MIDDLE\" BGCOLOR=\"#FFFFFF\">
<TD WIDTH=\"40\"><PRE>" <> hstart <> ToString[nsing] <> hend <> "</PRE></TD>
<TD WIDTH=\"200\"><PRE>" <> If[nsing==nsing0, "<B><FONT COLOR=\"#FF0000\">", ""] <> print2a[sings[[nsing]],4] <> If[nsing==nsing0, "</FONT></B>", ""] <> "</PRE></TD>
<TD WIDTH=\"200\"><PRE>" <> hstart <> print2p[weights[[nsing]],3] <> hend <> "</PRE></TD>" // p;
If[nsing==1,
"<TD ROWSPAN=" <> ToString[msing] <> " WIDTH=\"340\">
<IMG SRC=\"singsq" <> ToString[nind] <> ".gif?" <> rndm <> "\" WIDTH=300 HEIGHT=300 ALT=\"Singularities of quadratic " <>
indtxt <> " approximant\"></TD>" // p];
"</TR>" // p;
, {nsing, msing}];
"</TABLE>" // p;
"<FONT SIze=\"-2\"><A HREF=\"#top\"><IMG SRC=\"../../../../iconz/top.gif\" WIDTH=13 HEIGHT=9 BORDER=0 ALT=\"Top of Page\" ALIGN=\"MIDDLE\"></A> <A HREF=\"#top\">Top of the page</A> " // p;
"<BR>" // p;
If[newfigure,
sings0 = Union[{sing0}, {Conjugate[sing0]}];
xylist0 = {Re[#], Im[#]} & /@ sings0;
xylist1 = {Re[#], Im[#]} & /@ sings;
scale = 1.5;
pltq0 = MultipleListPlot[xylist0,
SymbolShape -> shape[1],
SymbolStyle -> {{Thickness[0.008 scale], color[0]}},
TextStyle -> {FontFamily -> "Times New Roman", FontSize -> 12,
FontSubstitutions -> {"Math1" -> "Symbol"}},
PlotRange -> xyrange{{-1, 1}, {-1, 1}},
AspectRatio -> 1,
Frame -> True,
ImageSize -> 300,
DisplayFunction -> Identity];
scale = 0.9;
pltq1 = MultipleListPlot[xylist1,
SymbolShape -> shape[0],
SymbolStyle -> {{Thickness[0.008 scale], color[1]}},
TextStyle -> {FontFamily -> "Times New Roman", FontSize -> 12,
FontSubstitutions -> {"Math1" -> "Symbol"}},
PlotRange -> xyrange{{-1, 1}, {-1, 1}},
AspectRatio -> 1,
Frame -> True,
ImageSize -> 300,
DisplayFunction -> Identity];
pltq = Show[pltq0, pltq1, DisplayFunction -> If[ifexport, Identity, $DisplayFunction]];
If[ifexport,
outfile = ToFileName[dir, "singsq" <>ToString[nind] <> ".gif"];
Export[outfile, pltq, "GIF", ImageSize -> 300, ImageResolution -> 72,
ConversionOptions -> {"AnimationDisplayTime" -> 1, "Loop" -> True,
"ColorReductionPalette" -> 16, "ColorReductionDither" -> False}]];
];
, {nind,mind}];
"<HR>" <> rulersing <> "<HR>" // p;
"<TABLE><TR ALIGN=\"CENTER\">\n
<TD ALIGN=\"CENTER\" VALIGN=\"CENTER\"><A HREF=\"index.htm#sings3\">Plot of singularities</A></TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"TOP\"><IMG SRC=\"../../../../blank.gif\" WIDTH=32 HEIGHT=43 ALT=\"Blank\"></TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"CENTER\"><A HREF=\"../index.htm\">" <>
"<IMG SRC=\"../../../../iconz/molec1.gif\" WIDTH=36 HEIGHT=41 BORDER=0 ALIGN=\"ABSMIDDLE\" ALT=\"Molecule - icon for Allen-data\">" <>
"List of examples</A></TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"TOP\"><IMG SRC=\"../../../../blank.gif\" WIDTH=32 HEIGHT=43 ALT=\"Blank\"></TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"CENTER\"><A HREF=\"../proglist.htm\"><I>Mathematica</I> programs</A></TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"TOP\"><IMG SRC=\"../../../../blank.gif\" WIDTH=32 HEIGHT=43 ALT=\"Blank\"></TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"CENTER\"><A HREF=\"../../index.htm\">Work in UMassD</A></TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"TOP\"><IMG SRC=\"../../../../blank.gif\" WIDTH=32 HEIGHT=43 ALT=\"Blank\">
</TD>\n
<TD ALIGN=\"CENTER\" VALIGN=\"CENTER\"><A HREF=\"../../../reports.htm\">" <>
"<IMG SRC=\"../../../waste.gif\" WIDTH=32 HEIGHT=43 BORDER=0 ALIGN=\"ABSMIDDLE\" ALT=\"Waste icon\">Unpublished reports</A></TD>\n
</TR></TABLE><HR>" // p;
"<P><FONT SIZE=\"-2\">Designed by <A HREF=\"../../../../index.htm\">A. Sergeev</A></FONT>.</P></BODY></HTML>" // p;
close[singsf];
of = oexampf;
printprogr[sect]; |