Test of the economolized rational approximants (ERAs) Some preliminaries ClearAll["Global`*"] Needs["Statistics`LinearRe+gression`"] Needs["Calculus`Pade`"] Now we define our test function (func[x]) and the value of the parameter (x0) at which this function will be evaluated. The series expansion (numberfunc) is then calculated up to some order (maxOrder). For comparison purposes, (straightTable) displays the convergence of the power series. The function era[i,n] is the [n/n] economized rational approximant evaluated at x0. The parameter (i) is the tunable ERA parameter. The function era[i,n] with i=0 is equivalent to the [n/n] Pade approximant, so we make a table of these Pades so as to compare with the results from the ERAs we calculate shortly. In[2622]:= func[x_] := Exp[x]; alpha = 8864/1000; x0 = 8; maxOrder = 18; numPrecision = 22; exactValue = N[func[x0], maxOrder]; era[i_, n_] := EconomizedRationalApproximation[func[x], {x, {-i, i}, n, n}]; Maclaurin[n_] := Normal[ Series[func[x], {x, 0, n}]] MaclaurinTable = Table[Maclaurin[n], {n, 0, maxOrder}] /. x -> x0; PadeTable = Table[era[0, n], {n, 0, maxOrder/2}] /. x -> x0; eraTable[alpha_] := Table[era[alpha, n], {n, 0, maxOrder/2}] /. x -> x0; Print["Maclaurin series: ", N[MaclaurinTable, 18] ]; Print["Error in Maclaurin series: ", N[ exactValue - MaclaurinTable, 18] ]; Print["Pade series: ", N[PadeTable, 18] ]; Print["Error in Pade series: ", N[exactValue - PadeTable, 18] ]; Print["ERA series: ", N[eraTable[alpha], 18] ]; Print["Error in ERA series: ", N[exactValue - eraTable[alpha], 18] ]; Now we calculate the sequence of ERAs, with the tunable parameter (i) being incremented in steps of (istep) in the outer Do[ ] loop. The inner While[ ] loop selects the best power law equation --- expressed as an exponential Exp[-k n] --- to fit the resulting table of values (notice that we don't try to fit the entire ERA sequence, only those values corresponding to the highest four orders). The function Regress[ ] [[2]] [[2]] calculates the regression for the chosen value of k. This While[ ] steps through a range of (k), starting with (kI), in steps of kStep. 1582, 12.3-17.1, 1273, 8856 In[2709]:= alphaI = 8856/1000; alphaStep = 1/1000; alphaF = 8856/1000; kI = 3.1; kstep = 1/10; kF = 3.1; fitOrder = 2; fitConstant = N[ 1 - 10^(-fitOrder)]; fitEquation[k_] := Fit[sequence, {1, Exp[-k n]}, n]; regression[k_] := Regress[sequence, {1, Exp[-k n]}, n][[2]][[2]]; ShowOutput[alpha_, k_] := Print["A(", fitOrder, ") ", N[alpha, 4], " ", N[k] , " F(n) = ", fitEquation[k], " ", MantissaExponent[N[1 - regression[k], 4]], " ", N[sequence , numPrecision] ] ; alpha = alphaI; While[ alpha <= alphaF, sequence = Take[eraTable[alpha], -4]; ListPlot[ sequence, PlotRange -> All, PlotJoined -> True]; k = kI; regTemp = 0; While[ k <= kF, If[ regression[k] > fitConstant, ShowOutput[alpha, k] ]; k += kstep; regTemp = regression[k]; ]; alpha += alphaStep ]; Print["Complete"] ; i = 8856/1000; n = 2; erat[i_, n_] := EconomizedRationalApproximation[numberfunc, {x, {-i, i}, n, n}] ; Do[a = NSolve[Denominator[erat[i, n]] == 0, x]; Print[a]; plota = Table[{Re[x /. a[[j]]], Im[x /. a[[j]]]}, {j, 1, Length[a]}]; Print["poles = ", plota]; b = NSolve[Numerator[erat[i, n]] == 0, x]; plotb = Table[{Re[x /. b[[j]]], Im[x /. b[[j]]]}, {j, 1, Length[b]}]; Print["zeroes = ", plotb]; Print[i]; alp = ListPlot[plota, DisplayFunction -> Identity, PlotRange -> {{0, 18}, {-14, 14}}]; blp = ListPlot[plotb, DisplayFunction -> Identity, PlotRange -> {{0, 18}, {-14, 14}}]; Show[alp, blp, DisplayFunction -> $DisplayFunction]; n += 1, {10}] x0 = 2; exactValue = Exp[8]; {x0, FortranForm[ exactValue - ERATable[[x0/2]]] } >> exptest2.math; Do[ {x0, FortranForm[ exactValue - ERATable[[x0/2]]] } >>> exptest2.math;; x0 += 2, {12}] Remove[g] g[y_] := N[ EconomizedRationalApproximation[numberfunc, {x, {-0, 0}, 4, 4}] /. x -> y, numprecision] x0 = -1.0; {x0, FortranForm[10^16(Exp[8] - g[x0])] } >> figone1.math; Do[ {x0, FortranForm[10^16(Exp[x0] - g[x0])] } >>> figone1.math;; x0 += .02, {110}] h[y_] := N[ EconomizedRationalApproximation[numberfunc, {x, {-1, 1}, 4, 4}] /. x -> y, numprecision] x0 = -1.03; {x0, FortranForm[ 10^16(Exp[8] - h[x0])] } >> figone2.math; Do[ {x0, FortranForm[ 10^16 (Exp[x0] - h[x0])] } >>> figone2.math;; x0 += .0, {110}] Print["PadeError ", errorTable]; plot1 = LogListPlot[errorTable, PlotRange -> {.000000001, 200000}, PlotJoined -> True, DisplayFunction -> Identity] ERATable = Table[N[era[10064/1000, n], numprecision], {n, 0, nmax/2 - 1}]; Print[nmax/2 - 1] Print["ERA ", N[ERATable, 20]]; errorTable = Abs[exactValue - ERATable]; Print["ERA Error: ", errorTable] plot2 = LogListPlot[errorTable, PlotRange -> {.000000001, 200000}, PlotJoined -> True, DisplayFunction -> Identity] Show[plot1, plot2, DisplayFunction -> $DisplayFunction] Print["PadeError ", errorTable]; plot1 = LogListPlot[errorTable, PlotRange -> {.000000001, 200000}, PlotJoined -> True, DisplayFunction -> Identity] ERATable = Table[N[era[10064/1000, n], numprecision], {n, 0, nmax/2 - 1}]; Print[nmax/2 - 1] Print["ERA ", N[ERATable, 20]]; errorTable = Abs[exactValue - ERATable]; Print["ERA Error: ", errorTable] plot2 = LogListPlot[errorTable, PlotRange -> {.000000001, 200000}, PlotJoined -> True, DisplayFunction -> Identity] Show[plot1, plot2, DisplayFunction -> $DisplayFunction]