# Plotting the energy function found by summation of the series

## Mathematica program

 ```(* Plotting the energy function found by summation of the series *) sect = "function"; entity = "figure"; (*Options*) mindlist = 10; (*number of approximants*) ncrop = 20; (*disregard coefficients with N > ncrop*) pltpnts = 999; If[iftest, pltpnts = 111]; f0 = func[[1]]; (* e0true = 0/10 f0; tr[e_, z_] := z e + (1 - z)e0true; *) tr[e_, z_] := e; cf = If[nm > ncrop, Take[func, ncrop], func]; 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 < 1 || nuse > nmax, Continue[]]; If[Abs[n1 - n2] > 1 || Abs[n2 - n3] > 1 || Abs[n3 - n1] > 1, Continue[]]; If[! (n1 <= n2 <= n3 <= n1 + 1), Continue[]]; indlist = Append[indlist, {ind, prior, nuse}], {n2, n1, n1 + 1}, {n3, n1, n1 + 1}], {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]] &]; tind = Transpose[indlist]; pind = {ToString /@ tind[[1]], tind[[3]]}; Print["Indexes: ", pind // TableForm]; (* Choosing range of z and E *) If[NumberQ[zmin] && NumberQ[zmax], ifzbounds = True; {zmin0, zmax0} = {zmin, zmax}, ifzbounds = False]; If[NumberQ[fmin] && NumberQ[fmax], iffbounds = True; {fmin0, fmax0} = {fmin, fmax}, iffbounds = False]; (* Approx. radius of convergence *) zr = 1/asave[name]; (* Approx. behavior of the function *) nf = 20; {zmin, zmax} = 0.85{-zr, -zr/3}; ztab = Table[zmin + (n - 1)/(nf - 1)(zmax - zmin), {n, nf}]; ftab = sumser[#, nm] & /@ ztab; zftab = Transpose[{ztab, ftab}]; ffitm = Fit[zftab, {1, z}, z]; {zmin, zmax} = 0.85{zr/3, zr}; ztab = Table[zmin + (n - 1)/(nf - 1)(zmax - zmin), {n, nf}]; ftab = sumser[#, nm] & /@ ztab; zftab = Transpose[{ztab, ftab}]; ffitp = Fit[zftab, {1, z}, z]; (* z - range *) slistr = Select[slist0, Head[#] === Real &]; slistrm = Select[slistr, # < 0 &]; slistrp = Select[slistr, # > 0 &]; slistc = Select[slist0, Head[#] === Complex &]; If[ifzbounds, {zmin, zmax} = {zmin0, zmax0}, zmin = Min[2.4 slistrm, 1.8 Re[slistc]]; If[Abs[zmin]<3 zr, zmin = 1.4 zmin]; If[zmin > 0, zmin = -5 zr]; zmax = Max[1.5 slistrp, 1.5 Re[slistc]]; If[Abs[zmax]<2 zr, zmax = 1.4 zmax]; If[zmax < 0, zmax = 5 zr]; ]; (* f - range *) If[iffbounds, {fmin, fmax} = {fmin0, fmax0}, {zmin0, zmax0} = {zmin, zmax}; {zmin, zmax} = {zmin0, 0}; ztab = Table[zmin + (n - 1)/(nf - 1)(zmax - zmin), {n, nf}]; ftab = (ffitm /. z -> #) & /@ ztab; fmax = Max[ftab]; {zmin, zmax} = {0, zmax0}; ztab = Table[zmin + (n - 1)/(nf - 1)(zmax - zmin), {n, nf}]; ftab = (ffitp /. z -> #) & /@ ztab; fmin = Min[ftab]; {zmin, zmax} = {zmin0, zmax0}; (* Streching f - range *) {fmin, fmax} = f0 + {1.4, 1.0} ({fmin, fmax} - f0); (* Including en + engap into f - range *) enex2 = enex + engap; fd = 0.07(fmax - fmin); If[NumberQ[engap] && enex2 + fd > fmax, fmax = enex2 + fd]; ]; If[iprint, Print[{{"Min","Max"},{zmin, zmax}, {fmin, fmax}} // TableForm]]; If[newfigure == True, scale = 1; dat = If[NumberQ[engap], {{{1, enex}}, {{1, enex2}}}, {{{1, enex}}}]; pltex = MultipleListPlot[Sequence @@ dat, SymbolShape -> {shape[0], shape[1]}, SymbolStyle -> {{Thickness[0.004 scale], RGBColor[1, 0, 0]}}, DisplayFunction -> Identity]; grat = 1/GoldenRatio // N; pltseq = {}; Clear[pol, c, funs]; (*Cycle over indexes*) Do[ {{n1, n2, n3}, prior, nc} = indlist[[nindlist]]; plts = Plot[tr[sumser[z, nc], z], {z, zmin, zmax}, PlotRange -> {{zmin, zmax}, {fmin, fmax}}, PlotStyle -> {Thickness[0.005], RGBColor[1, 0, 0]}, AxesOrigin -> {0, f0}, PlotPoints -> pltpnts, DisplayFunction -> Identity]; padefrac = padesum[z, nc]; padefun[t_] := ( t1 = SetPrecision[t, 64]; ft = padefrac /. z -> t1; If[ft > fmin && ft < fmax, ft, "near Pole"]); Off[Plot::"plnr"]; pltp = Plot[tr[padefun[z], z], {z, zmin, zmax}, PlotRange -> {{zmin, zmax}, {fmin, fmax}}, PlotStyle -> {Thickness[0.005], RGBColor[0, 0, 1]}, AxesOrigin -> {0, f0}, PlotPoints -> pltpnts, Compiled -> False, DisplayFunction -> Identity]; 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]]; polfz = pattern /. s1 /. f[z] -> f; pola = pol[2] /. s1; dpola = D[pola, z]; Clear[root]; dim = 10.^-32; s = Solve[polfz == 0, f]; dselect = 0.1(zmax - zmin)/pltpnts; root[t_, n_] := (t1 = SetPrecision[t, 64]; fnt = (f /. s[[n]] /. z -> t1); {refnt, imfnt} = {Re[fnt], Im[fnt]}; If[Abs[imfnt/fnt] < dim, If[Abs[(pola/dpola) /. z -> t1] > dselect, refnt, "almost Infinity"], "complex-valued"]); reroot[t_, n_] := (t1 = SetPrecision[t, 64]; fnt = (f /. s[[n]] /. z -> t1); {refnt, imfnt} = {Re[fnt], Im[fnt]}; If[Abs[imfnt/fnt] > dim, If[Abs[(pola/dpola) /. z -> t1] > dselect, refnt, "almost Infinity"], "real-valued"]); pltq = Plot[{tr[root[z1, 1], z1], tr[root[z1, 2], z1]}, {z1, zmin, zmax}, PlotRange -> {{zmin, zmax}, {fmin, fmax}}, PlotPoints -> pltpnts, PlotStyle -> {{Thickness[0.008], RGBColor[0, 1, 0]}}, AxesOrigin -> {0, f0}, DisplayFunction -> Identity ]; pltqre = Plot[{tr[reroot[z1, 1], z1], tr[reroot[z1, 2], z1]}, {z1, zmin, zmax}, PlotRange -> {{zmin, zmax}, {fmin, fmax}}, PlotPoints -> pltpnts, PlotStyle -> {{Thickness[0.008], RGBColor[0, 1, 0], Dashing[{0.005, 0.015}]}}, AxesOrigin -> {0, f0}, DisplayFunction -> Identity ]; On[Plot::"plnr"]; pltg = Show[pltqre, pltq, pltp, plts, pltex, PlotRange -> {{zmin, zmax}, {fmin, fmax}}, AxesOrigin -> {0, f0}, TextStyle -> {FontFamily -> "Times New Roman", FontSize -> 12, FontSubstitutions -> {"Math1" -> "Symbol"}}, AspectRatio -> grat, Frame -> True, FrameLabel -> {"\!\(\* StyleBox[\"z\",\nFontFamily->\"Times\",\nFontSize->16,\nFontSlant->\"Italic\"]\)", "\!\(\* StyleBox[\"E\",\nFontFamily->\"Times\",\nFontSize->16,\nFontSlant->\"Italic\"]\)"}, ImageSize -> 400, DisplayFunction -> Identity]; (*Show together*) wr1 = 0.04; wr2 = 0.06; wr3 = 0.087; wd = 0.006; {ybot, ytop} = {0.122, 0.987} + 0.0{1, -1}; pltr = Rectangle[{0, 0} // Scaled, {1 - wr1 - wr2 - wr3 - 3 wd, 1} // Scaled, pltg]; t1 = Text["Sum", {1/2, 1/2}, TextStyle -> {FontFamily -> "Times New Roman", FontSize -> 12, FontWeight -> "Bold", FontColor -> RGBColor[1, 0, 0]}]; t2 = Text["Pade", {1/2, 1/2}, TextStyle -> {FontFamily -> "Times New Roman", FontSize -> 12, FontWeight -> "Bold", FontColor -> RGBColor[0, 0, 1]}]; t3 = Text["Quadr.", {1/2, 1/2}, TextStyle -> {FontFamily -> "Times New Roman", FontSize -> 12, FontWeight -> "Bold", FontColor -> RGBColor[0, 1, 0]}]; g = Graphics[#, Background -> RGBColor[0.9, 0.9, 1]] & /@ {t1, t2, t3}; mindlist1 = mindlist + 1; y1 = ybot + (ytop - ybot)(mindlist1 - 1)/mindlist1; y2 = ybot + (ytop - ybot)(mindlist1)/mindlist1; x2 = 1; x1 = x2 - wr3; r3 = Rectangle[{x1, y1} // Scaled, {x2, y2} // Scaled, g[[3]]]; x2 = x1 - wd; x1 = x2 - wr2; r2 = Rectangle[{x1, y1} // Scaled, {x2, y2} // Scaled, g[[2]]]; x2 = x1 - wd; x1 = x2 - wr1; r1 = Rectangle[{x1, y1} // Scaled, {x2, y2} // Scaled, g[[1]]]; pltindh = {r1, r2, r3}; pltind = Table[{{n1, n2, n3}, prior, nc} = indlist[[n]]; np1 = Floor[(nc - 1)/2]; np2 = nc - 1 - np1; txt = If[n <= mind, {ToString[nc], "[" <> ToString[np1] <> "/" <> ToString[np2] <> "]", "[" <> ToString[n1] <> ", " <> ToString[n2] <> ", " <> ToString[n3] <> "]"}, {"", "", ""}]; t = Text[#, {1/2, 1/2}, TextStyle -> {FontFamily -> "Times New Roman", FontSize -> 14, FontWeight -> If[n == nindlist, "Bold", "Plain"], FontColor -> RGBColor[0, 0, 0.5]}] & /@ txt; col = If[n > mind, RGBColor[1, 1, 1], If[n == nindlist, RGBColor[1, 1, 0], RGBColor[3/4, 3/4, 3/4]]]; g = Graphics[#, Background -> col] & /@ t; y1 = ybot + (ytop - ybot)(mindlist1 - (n + 1))/mindlist1; y2 = ybot + (ytop - ybot)(mindlist1 - (n + 1) + 1)/mindlist1; x2 = 1; x1 = x2 - wr3; r3 = Rectangle[{x1, y1} // Scaled, {x2, y2} // Scaled, g[[3]]]; x2 = x1 - wd; x1 = x2 - wr2; r2 = Rectangle[{x1, y1} // Scaled, {x2, y2} // Scaled, g[[2]]]; x2 = x1 - wd; x1 = x2 - wr1; r1 = Rectangle[{x1, y1} // Scaled, {x2, y2} // Scaled, g[[1]]]; {r1, r2, r3}, {n, mindlist}]; plt = Show[{pltr, pltindh, pltind} // Graphics, AspectRatio -> (1 - wr1 - wr2 - wr3 - 3 wd)grat, ImageSize -> 660, DisplayFunction -> If[ifexport, Identity, \$DisplayFunction]]; If[ifexport, outfile = ToFileName[dir, sect <>ToString[nindlist] <> ".gif"]; Export[outfile, plt, "GIF", ImageSize -> 660, ImageResolution -> 72, ConversionOptions -> {"AnimationDisplayTime" -> 1, "Loop" -> True, "ColorReductionPalette" -> 32, "ColorReductionDither" -> False}]; pltseq = Append[pltseq, plt] ]; , {nindlist, mind}]; If[ifexport, outfile = ToFileName[dir, sect <> ".gif"]; Export[outfile, pltseq, "GIF", ImageSize -> 660, ImageResolution -> 72, ConversionOptions -> {"AnimationDisplayTime" -> 1, "Loop" -> True, "ColorReductionPalette" -> 32, "ColorReductionDither" -> False}]; ]; ]; " sect <> "\">
" // p; "" // p; "" // p; "" <> "
The function E(z) found by summation of its power series.
Dashed line indicates that the approximant is complex valued.
Red dot marks exact physical energy at z = 1." // p; If[NumberQ[engap], "
Red circle marks the lowest excited energy level at z = 1." // p ]; "
To view results of summation of a specific number of terms of the series, click on the right bar.
" <> " sect <> ".gif?" <> rndm <> "\" WIDTH=660 HEIGHT=324 BORDER=0 ALT=\"Partial sums, " <> "Pade and quadratic approximants\" USEMAP=\"#indbar1\">
" // p; x0 = 527; x1 = 658; ytop = 28; ybot = 283; Do[y1 = ytop + n/mindlist (ybot-ytop); y0 = ytop + (n-1)/mindlist (ybot-ytop); {x0p, x1p, y0p, y1p} = ToString/@Round/@{x0, x1, y0, y1}; ref = sect <> ToString[n] <> ".gif?" <> rndm; {{n1, n2, n3}, prior, nc} = indlist[[n]]; app = "[" <> ToString[n1] <> ", " <> ToString[n2] <> ", " <> ToString[n3] <> "]"; " x0p <> ","<> y0p <> ","<> x1p <> ","<> y1p <> "\" href=\"" <> ref <> "\" onclick=\"vappr=window.open('" <> ref <> "','Summation" <> (*ToString[n]*)"" <> "','toolbar=no,status=no,scrollbars=no, location=no,menubar=no,directories=no,width=670,height=334'); vappr.focus(); return false\">" // p; , {n, mind}]; "" // p; printprogr[sect]; printnav;```

 Examples of MP seriesMathematica programsWork in UMassDUnpublished reports

Designed by A. Sergeev.