Message Boards Message Boards

[GIF] Plotting the Contours of Deformed Hyperspheres

Posted 7 years ago

Considered as a series reversion ( Cf. Mathworld, A&S ) of implicit equation $$E = \frac{1}{2}\Psi^2 + \sum_{n=3}^{\infty} f_n(\phi,\theta,\ldots)\;\Psi^n,$$ function A276738 determines the radius $\Psi$ of a hypersurface that limits in small $E$ to the shape of a perfect hypersphere, a circle, a sphere, etc...

Examples in 2,3 dimensions lend themselves well to depiction. Let's see the Henon Heiles Potential, and an octahedral Energy Surface.

Basic Functions

RExp[n_] := Expand[b Plus[R[0], Total[b^# R[#] & /@ Range[n]]]]

RCalc[n_] := 
With[{basis = Subtract[Tally[Join[Range[n + 2], #]][[All, 2]],
Table[1, {n + 2}]] & /@ IntegerPartitions[n + 2][[3 ;; -1]]},
Total@ReplaceAll[Times[-2, Multinomial @@ #, v[Total[#]],
Times @@ Power[RSet[# - 1] & /@ Range[n + 2], #]] & /@ basis,
{Q^2 -> 1, v[2] -> 1/4}]]

AbsoluteTiming[RSet[0] = 1;   Set[RSet[#], Expand@RCalc[#]] & /@ Range[20];][[1]]

Surface[Nexp_, rep_] :=   RExp[Nexp] /. R -> RSet /. 
v[n_] :>  Function[{a}, 
Total[v[#, a - #] Q^# P^(a - #) & /@ Range[0, a]]][n] /. 
rep /. {P -> Sin[\[Phi]], Q -> Cos[\[Phi]]};

HyperSurface[Nexp_, rep_] :=   RExp[Nexp] /. R -> RSet /. v[n_] :>    Function[{a}, 
Total[Flatten@          Table[v[a - i - j, i, j] X^(a - i - j) Y^i Z^j, {i, 0,  a}, {j, 0, a - i}]]][n] /. 
rep /. {X -> Sin[\[Theta]] Cos[\[Phi]], 
Y -> Sin[\[Theta]] Sin[\[Phi]], Z -> Cos[\[Theta]]};

Print Pictures

ViewV = {1, 2, 2};

TrigWords = 
Flatten[Expand[      Surface[20, {v[2, 1] -> 3 \[Epsilon], v[0, 3] -> -\[Epsilon], 
v[_, _] -> 0}]] /. Plus -> List /. Times -> List];

TrigLines20 = 
With[{surf =     Surface[20, {v[2, 1] -> 1/2, v[0, 3] -> -1/3/2, v[_, _] -> 0}]},
(surf /. b -> (#/10)) & /@ Range[5]];

g1 = Show[
Graphics[{Lighter@Gray,  Text[#, RandomReal[{-1/2, 1}, 2]] & /@ TrigWords}, 
PlotRange -> {{-1/2, 1}, {-1/2, 1}}], Graphics[{Thick,
{Dashed,  Line[Part[RotationMatrix[{{0, 0, 1}, ViewV}].#, 
1 ;; 2] & /@ {{0, 0, -1/4}, {0, 0, 5/4}}]},
Thickness[.005],  MapThread[
Line[Function[{a},     N[Part[RotationMatrix[{{0, 0, 1}, 
ViewV}].{#1 Cos[\[Phi]], #1 Sin[\[Phi]], #2/
10} /. {\[Phi] -> a}, 1 ;; 2]]] /@ (Range[0, 200]/
200 2 Pi)] &,
{TrigLines20, Range[5]}]}], ImageSize -> 800];

OctSurf6 =   HyperSurface[
20, {v[4, 0, 0] -> \[Epsilon] 2/5, v[0, 4, 0] -> \[Epsilon] 2/5, 
v[0, 0, 4] -> \[Epsilon] 2/5,
v[2, 2, 0] -> \[Epsilon] 6/5, v[2, 0, 2] -> \[Epsilon] 6/5, 
v[0, 2, 2] -> \[Epsilon] 6/5, 
v[_, _, _] -> 0} /. {\[Epsilon] -> -1/2}] /. b -> 6/10;

OctSurf3 =   HyperSurface[
20, {v[4, 0, 0] -> \[Epsilon] 2/5, v[0, 4, 0] -> \[Epsilon] 2/5, 
v[0, 0, 4] -> \[Epsilon] 2/5,
v[2, 2, 0] -> \[Epsilon] 6/5, v[2, 0, 2] -> \[Epsilon] 6/5, 
v[0, 2, 2] -> \[Epsilon] 6/5, 
v[_, _, _] -> 0} /. {\[Epsilon] -> -1/2}] /. b -> 3/10;

OctWords =   Flatten[Expand[
HyperSurface[  10, {v[4, 0, 0] -> \[Epsilon] 2, v[0, 4, 0] -> \[Epsilon] 2, 
v[0, 0, 4] -> \[Epsilon] 2 ,
v[2, 2, 0] -> \[Epsilon] 6 , v[2, 0, 2] -> \[Epsilon] 6 , 
v[0, 2, 2] -> \[Epsilon] 6 , v[_, _, _] -> 0}]] /. 
Plus -> List /. Times -> List];

g2 = Show[   Graphics[{Lighter@Gray, 
Text[#, RandomReal[{-2, 2}, 2]] & /@ OctWords}, 
PlotRange -> 2 {{-1, 1}, {-1, 1}}],
Graphics[{     Thick, {Dashed, 
Line[Part[RotationMatrix[{{0, 0, 1}, ViewV}].#, 1 ;; 2] & /@ {{0, 0, 2}, {0, 0, -2}}]},
Thickness[.005],
Line /@       Outer[N[Part[
RotationMatrix[{{0, 0, 1}, 
ViewV}].{Sin[\[Theta]] Cos[\[Phi]], 
Sin[\[Theta]] Sin[\[Phi]], 
Cos[\[Theta]]} OctSurf3 /. {\[Theta] -> #1, \[Phi] ->  #2}, 1 ;; 2]] &, 
Range[1, 9] Pi/10, Range[0, 100] 2 Pi/100, 1],
Line /@       Outer[N[Part[
RotationMatrix[{{0, 0, 1}, 
ViewV}].{Sin[\[Theta]] Cos[\[Phi]], 
Sin[\[Theta]] Sin[\[Phi]], 
Cos[\[Theta]]} OctSurf6 /. {\[Theta] -> #1, \[Phi] ->   #2}, 1 ;; 2]] &, 
Range[1, 9] Pi/10, Range[0, 200] 2 Pi/200, 1]
}], ImageSize -> 800];

Henon Heiles

Henon Heiles

Octahedral

enter image description here

What's more curious than these graphics? The new conjecture: that A276738 can be defined in terms of A028338.

POSTED BY: Brad Klee
13 Replies

To be clear... around 1:00 a.m. 3/7/2017 I did forget to type out the definition of RExp, and later edited it in today.

Temper, ha! My current preference is distemperment. To each his own.

To sum things up, one more point :

  • If the built-in functions don't meet user specifications, it's easy for users to code their own improvements using the basic functionality of Mathematica.

Thanks again for your comments, I will keep working on improving the definitions for these algorithms.

Brad

POSTED BY: Brad Klee

Hi Danny,

Sorry to say this, but it seems like you are not reading the posts at all. A few more points:

  • You say that RExp is not defined. It is defined in two previous posts. I copied the definition again into the post from last night. Here it is a fourth time:

    RExp[n_] := Expand[y Plus[R[0], Total[y^# R[#] & /@ Range[n]]]]
    
  • You suggest that I look at mathworld, but I already have referenced that website in two previous posts.

As last night's post says, mathworld does provide a solution for $m=1$. Furthermore, this solution does not apply to any other $m$ as the $a_1$ coefficient appears in the denominator of all $A_i$, so cannot equal to $0$. I checked Morse and Feshbach and did not see anything beyond $m=1$. But Eq. 12 in the Mathworld article is what I call the Coefficient Grading Condition. Notice that the number of integer solutions to this equation exactly equals the partition numbers. I suspect that the generalization of Eq. 11 to arbtirary $m$ is exactly equivalent to the Multinomial method I present here.

  • Your latest code is a nice simplification, but not really different than mine.

Without loss of generality we can assume $y>0$, so I added Refine. Also if you plan to compare using SameQ, you need to use Expand rather than PowerExpand . I'm now using:

mSeries2[n_, m_] :=  x^m + Array[c, n, m + 1].x^Range[m + 1, m + n] + O[x]^(m + n + 1)
mExpand2[n_, m_] :=  Refine[Expand[Normal[InverseSeries[mSeries2[n, m], z]] /. z -> y^m], y > 0]

All the time dependence is tied up in InverseSeries, so we see the same behavior as above. For $n=12$ the calculation takes about $30(s)$, which suggests that $n=13$ takes over $100(s)$, $n=14$ over $400(s)$.

  • You say: <<As for having a "general form" applicable to the case of independent symbolic coefficients, it would not surprise me if a direct method applies, and is substantially more efficient. >>

As the above time graph shows, there are not one but two direct methods, both substantially more efficient ! The multinomial method is easy to understand, but I think the second mysterious number theory method will take some more effort to understand and prove.

  • Concession: I mispoke about the complexity class. Computation of expansion term $n$ is $\mathcal{O}(p(n))$. Computation of the entire series should then be $\mathcal{O}(\sum_{i}^{n} p(i))$.

  • In general, the dilligence of referees is dissapointing me, and I'm getting tired of what nowadays passes for acceptable criticism.

  • After a few minutes I get tired of wating on the following command:

    In[78]:= AbsoluteTiming[mExpand2[15, 2]]
    Out[78]= $Aborted
    
  • Again, Compare with Multinomial algorithm:

    In[80]:= AbsoluteTiming[MultinomialExpand[15, 2]][[1]]
    Out[80]= 0.436722
    
  • And then ( in the fanfiction voice ):

    (* LORD SO-AND-SO: << Mein zauberstab, schnell !  I'm feeling like a magical supremacist again. >> *)
    AvadaKedavra = AbsoluteTiming[MagicExpand[15, 2]][[1]]
    Out[]=0.055286 
    

AvadaKedavra

POSTED BY: Brad Klee

Right you are, I cut-and-pasted and apparently missed the first line. But expecting anyone to read carefully through all prior posts strikes me as a reach: I've spent hours writing responses for forums, with no expectation that more than 20% might be read by anyone. Then again, I've been doing this for a while and have tempered my expectations.

As for the m>1 case, it can be reduced by inverting the product of x times the mth root of a unit (a series of the form 1+x*something). Whether that gives rise to the form you have I do not know offhand.

POSTED BY: Daniel Lichtblau

Hi, Danny, thanks for taking some time to respond. I think the code is already out there to satisfy all of your demands, but admittedly I have not been explaining the calculations in any detail. Let's aim for closer to article quality with this post, and see if we can convince ourselves that Mathematica function InverseSeries can be improved on a wide range of inputs. We have not two, but three distinct methods to compare!

Statement of the Problem & Variable Parameters

We start with a power series expansion $$ z = x^m+\sum_{i>m}^{\infty}c_i \; x^i$$ with $m$ defined as the leading power. Without loss of generality we have already assumed that the coefficient of the leading term $c_m$ can be set to $1$ by change of the $x$ variable. We seek a solution $x=f(y(z))$, which expands in whole powers of $y(z)$,

$$ x = y + \sum_{i>1}^{n+1} f_{i}(\mathbf{c}) y^i, $$

with $n$ defined as truncation parameter and $f_{i}(\mathbf{c})$ is in general a complicated, but well structured function of the expansion coefficients $\mathbf{c}= ( c_i:i\in1,2,\ldots )$. It's worthwhile to mention up front that the $f_{i}(\mathbf{c})$ for case $m=1$ have already been explored in numerous outlets including Mathworld, Abramowitz and Stegun, and the OEIS. From the OEIS entry: "The sequence of row lengths is A000041(n) (partition numbers)", which gives us an expectation regarding the number of terms in each $f_i(\mathbf{c})$. This same row-length feature is shared by the triangle for $m=2$, which I added to the OEIS while working on physics problems, including planetary motion, A276738. This is an important feature at any $m$, which we will return to shortly. All functions will take parameters $m$ and $n$ as input and return as output the expansion $x(y(z))$, which inverts the $z(x)$ expansion. Our complexity expectations are as follows:

  • Changing $m$ does not change complexity class
  • Changing $n$ explores complexity class $\mathcal{O}(p(n))$, where the $p(n)$ are the partition numbers, listed in A000041.

Using the built-in Function

Let's take the code from above MmaExpand, and break it apart into two functions

MmaSeries[n_, m_] := Series[f[x], {x, 0, n + m}] /. 
Join[{f[0] ->  0, (D[f[x], {x, m}] /. x -> 0) -> (m!)}, (D[f[x], {x, #}] /. x -> 0) -> 0 & /@ 
Range[1, m - 1], ((D[f[x], {x, #}] /. x -> 0) -> #! c[#]) & /@ 
Range[m + 1, n + m]]

MmaExpand[n_, m_] := Expand[Refine[Normal[InverseSeries[MmaSeries[n, m]]] /. x -> y^m,  y > 0]]

The first function generates a series and the second function--It's magic to me, slower magic--generates the inverse expansion. Additionally we utilize the small $x$ limit where $z \approx x^m$ to prove that the solution is an expansion in integer powers of $y=z^{1/m}$. This is easy to see once you have studied the multinomial method in the next section. For now we need to run a few critical tests:

AbsoluteTiming[
 RowLengthTest = MatrixForm[Function[{m},
     Length /@ (Coefficient[MmaExpand[8, m], y, #] /. Plus -> List & /@
         Range[3, 8])] /@ Range[5]
   ]]

AbsoluteTiming[ Outer[Coefficient[  Expand[Normal@MmaSeries[4, #1] /. x -> MmaExpand[4, #1]], 
     y, #2] &,  Range[5], Range[5] ] // MatrixForm ]

out1

Both are maps where index $m=1,2,3,4,5...$ corresponds to a row. First, we see the number of terms is the partition numbers, as expected. Next in the identity matrix, we see up to power $y^5$ substitution of the inverse series into the series cancels all terms other than $y^m=z$, as required. Let's take for granted that the only problem is time. We obtain timing statistics by mapping over $m$ and $n$

t1=Outer[AbsoluteTiming[MmaExpand[#2, #1]][[1]] &, Range[5], Range[10]]
Out[] = {
{0.00795752, 0.0104603, 0.00681345, 0.00671655, 0.0118878, 0.022436,  0.0632026, 0.159615, 0.443813, 1.42085}, 
{0.00181332, 0.00326319,  0.00738307, 0.0140106, 0.0274177, 0.0521328, 0.118843, 0.26841, 0.684693, 1.9761}, 
{0.00188818, 0.00317625, 0.00734353, 0.013728,   0.0268191, 0.051075, 0.116821, 0.266313, 0.681015,  1.92303}, 
{0.00231411, 0.00382979, 0.0077761, 0.0146602, 0.0272598,  0.0552915, 0.119205, 0.270802, 0.690533, 1.93577},
{0.00201708,  0.00370995, 0.00782772, 0.014608, 0.0296672, 0.0543108, 0.119673, 0.269981, 0.702621, 1.92741}}

We compare these times with other methods.

Multinomial Method

as above:

RExp[n_] := Expand[y Plus[R[0], Total[y^# R[#] & /@ Range[n]]]]
RCalc[0, _] = 1;
RCalc[n_, m_] := With[{basis = Subtract[Tally[Join[Range[n + m], #]][[All, 2]], 
Table[1, {n + m}]] & /@   Select[IntegerPartitions[n + m],    Length[#] > m - 1 &][[2 ;; -1]]}, 
Total@ReplaceAll[    Times[-1/m, Multinomial @@ #, c[Total[#]], 
Times @@ Power[RSet[# - 1] & /@ Range[n + m], #]] & /@     basis, {c[m] -> 1}]]
MultinomialExpand[n_, m_] := (Clear@RSet;   Set[RSet[#], Expand@RCalc[#, m]] & /@ Range[0, n];
Expand[RExp[n] /. R[o_] :> RCalc[o, m]])

The code is a conceptually simple, straightforward realization of the solution method discussed for case $m=2$ in Plane Pendulum and Beyond (rejected by American Journal of Physics), generalized to any $m$. This method is easy to prove, especially for someone versed in ring theory. The algorithm takes advantage of a graded coefficient-algebra associated with partitions of the integer numbers. Rather than writing out a full proof, we can test equivalence:

SameQ[MultinomialExpand[10, #] & /@ Range[5],  MmaExpand[10, #] & /@ Range[5]]
Out[]= True

And now we can run the timing statistics up to $n=18$

t2 = Outer[AbsoluteTiming[MultinomialExpand[#2, #1]][[1]] &, Range[5], Range[18]]
Out[]={{0.000350166, 0.000335676, 0.0006466, 0.00102544, 0.00187701, 
  0.00335827, 0.0050605, 0.00829471, 0.0137534, 0.0213327, 0.0338478, 
  0.0526303, 0.0818673, 0.127584, 0.197492, 0.303197, 0.461779, 
  0.705808}, {0.00059166, 0.000385183, 0.000731424, 0.00141787, 
  0.00256708, 0.00438643, 0.0072424, 0.0121079, 0.0196304, 0.0317456, 
  0.0498247, 0.0797138, 0.123864, 0.19419, 0.297837, 0.463815, 
  0.703707, 1.08945}, {0.000783647, 0.000391522, 0.000853077, 
  0.00165484, 0.00296434, 0.00541067, 0.00877981, 0.0160512, 
  0.0248739, 0.0400204, 0.0637266, 0.101257, 0.158744, 0.247118, 
  0.392642, 0.600964, 0.915037, 1.42288}, {0.000755876, 0.000407823, 
  0.000880547, 0.00177226, 0.0032173, 0.00567963, 0.00968088, 0.01652,
   0.0290402, 0.0468616, 0.0753077, 0.120092, 0.191354, 0.298932, 
  0.466456, 0.734795, 1.09722, 1.70526}, {0.000888697, 0.000464272, 
  0.000976541, 0.00196093, 0.00343374, 0.00653271, 0.0109065, 
  0.0190267, 0.0313326, 0.0523583, 0.0833902, 0.127929, 0.202708, 
  0.327367, 0.512287, 0.798683, 1.23594, 1.91088}}

Still we are under two seconds.

Magic Number Theory Method

as above:

SF[0, 0, m_] = 1;SF[n_, k_, m_] := SF[n, k, m] = 
If[Or[k > n, k < 0], 0,    If[And[n == 0, k == 0], 1, 
SF[n - 1, k - 1, m] + (m*n - 1)*SF[n - 1, k, m]]]
GeneralPolys[n_, m_] :=  MapThread[  Dot, {Table[SF[i, j, m], {i, 0, n - 1}, {j, 0, i}],    
Table[x^j, {i, 0, n - 1}, {j, 0, i}]}]
GeneralT[n_, m_] :=  With[{ps = GeneralPolys[n, m]},   
Table[(-1)^j ps[[j]] /. {x -> i + 2}, {i, 1, n}, {j, 1, i}]]

Basis[n_, m_] :=  Total /@ Map[
Times[m^(-Length[#]) 1/Times @@ (Tally[#][[All, 2]]!), 
Times @@ (# /. x_Integer :> c[x + m])] &, 
Function[{a}, Select[IntegerPartitions[n], Length[#] == a &]] /@ 
Range[n], {2}]

MagicExpand[n_, m_] := Expand[y +   Dot[MapThread[Dot, {GeneralT[n, m], Basis[#, m] & /@ Range[n]}], 
y^Range[2, n + 1]]];

Exploring the expansion coefficients in the OEIS partition triangles A111785 and A276738, I conjectured that it would be possible to store only a regular triangle, i.e., row length equals row index. For example, see: A283247, and triangles for $m=1,2,3...12$ are given in a previous post above. For the final formulation I used Peter Luschny's definitions for the Stirling-Frobenius Cycle numbers . It's conceivable that the above algorithm is easy for a pure mathematician to prove, and also that the above algorithm is already known to someone. I don't know. I found this method by a smart empirical investigation. At this time I can only say that it works, relative to other proven methods, to many orders of magnitude. Here we can check $15$ quickly

SameQ[MagicExpand[15, #] & /@ Range[15],  MultinomialExpand[15, #] & /@ Range[15]]
Out[72]= True

This is a promising result. Although I'm not trained as a mathematician, and mathematicians don't seem to like me very well lately, I'm interested to work on a proof starting from the known algorithm using multinomial coefficients. For now, let's just do the timing test, up to $n=30$

t3 = Outer[AbsoluteTiming[MagicExpand[#2, #1]][[1]] &, Range[5], Range[30]]

Out[75]= {{0.000517401, 0.000730217, 0.00129441, 0.00241826, 
  0.00313761, 0.00480814, 0.00677571, 0.00235547, 0.00322757, 
  0.00458355, 0.00661482, 0.00910492, 0.0127276, 0.0178956, 0.0245657,
   0.03328, 0.045086, 0.0608486, 0.0816343, 0.111034, 0.145875, 
  0.192592, 0.252839, 0.331296, 0.431803, 0.554355, 0.709992, 
  0.907644, 1.15345, 1.45812}, {0.000193195, 0.000157575, 0.000248437,
   0.000478158, 0.000638147, 0.000984389, 0.00152262, 0.00219397, 
  0.00318681, 0.00455035, 0.00642736, 0.00907141, 0.0128466, 
  0.0179599, 0.0244666, 0.0334372, 0.0457468, 0.0618206, 0.0836055, 
  0.115272, 0.151878, 0.196921, 0.259114, 0.338644, 0.436268, 
  0.564819, 0.728962, 0.920223, 1.17813, 1.47938}, {0.000208892, 
  0.000162405, 0.000285868, 0.000402993, 0.000645996, 0.00100673, 
  0.00151145, 0.00223261, 0.00324719, 0.00460649, 0.00675066, 
  0.00924921, 0.0127491, 0.0177939, 0.0250812, 0.0334789, 0.0458992, 
  0.0638416, 0.0865447, 0.112189, 0.149649, 0.198407, 0.25906, 
  0.34008, 0.436906, 0.56552, 0.727684, 0.928567, 1.18617, 
  1.52182}, {0.000167536, 0.000164216, 0.000273794, 0.000405408, 
  0.000670145, 0.000981371, 0.00175083, 0.00223442, 0.00321972, 
  0.00472271, 0.00701268, 0.0105626, 0.0139251, 0.0198001, 0.0275424, 
  0.0371674, 0.0510624, 0.0657576, 0.0853008, 0.120408, 0.156297, 
  0.213033, 0.280924, 0.365606, 0.467034, 0.589213, 0.780763, 1.02078,
   1.23597, 1.59574}, {0.000196817, 0.000165725, 0.000252361, 
  0.000439821, 0.000629393, 0.000999784, 0.00157846, 0.0022308, 
  0.00322183, 0.00483742, 0.00662297, 0.00945479, 0.0128339, 
  0.0181477, 0.0254782, 0.034367, 0.0472682, 0.0632687, 0.0855504, 
  0.113633, 0.152079, 0.201328, 0.268378, 0.346395, 0.478177, 
  0.617492, 0.742203, 0.935041, 1.20908, 1.52182}}

again, under two seconds max time.

Graphing Time Results

Show[
 ListLinePlot[t1,  PlotStyle -> (Blend[{Red, Orange}, #/5] & /@ Range[5]), 
  PlotRange -> All], ListLinePlot[t2, PlotStyle -> (Blend[{Blue, Green}, #/5] & /@ Range[5]), 
  PlotRange -> All],  ListLinePlot[t3, PlotStyle -> Purple, PlotRange -> All], 
 PlotRange -> All, AxesLabel -> {Style["n", 20], Style["t(s)", 20]}, ImageSize -> 800 ]

Timing Test

This is the same graph as above, but with extra explanation, maybe it is easier to interpret. As predicted, variation with $m$ is insignificant compared to the blow up with $n$, associated with partition complexity class. Compare with a plot of the partition numbers:

ListLinePlot[Length@IntegerPartitions[#] & /@ Range[20], ImageSize -> 800]

Partition Numbers

The time graph also shows that in the partition complexity class, regardless of $m$, the magic method vastly outperforms other methods. The easy multinomial method only outperforms the method using built-in inverse function. Setting a time limit of $2s$, we can only compute approximately $10$, $18$, $30$ orders of magnitude, as indexed by $n$.

Again, the most interesting feature of the graph is the stability of the magic method with regard to variation of $m$. If you look closely at the code, tasks seperate into an $\mathcal{O}(n^2)$ listing of coefficients, and an $\mathcal{O}(p(n))$ expansion of functions $g_{i,j}(\mathbf{c})$, called a "basis". Again, if the basis function can be optimized by recursive programming, then the timing can be further reduced.

Danny, does this make more sense? If not maybe we can get a third opinion? Seems like Michael Trott would probably also have some idea about this sort of calculation. Ultimately it would also be interesting to hear from one of you inside the company: How does the Mma algorithm work for inverting a series with symbolic coefficients?

That's all for now,

Brad

POSTED BY: Brad Klee

My first remark is that the code for the multinomial method does not work for me. I think there may be a definition missing, for RExp.

If I follow correctly, your basic version might be shortened as below.

mSeries2[n_, m_] := 
 x^m + Array[c, n, m + 1].x^Range[m + 1, m + n] + O[x]^(m + n + 1)
mExpand2[n_, m_] := 
 PowerExpand[Normal[InverseSeries[mSeries2[n, m]]] /. x -> y^m]

Example:

In[30]:= mSeries2[5, 2]

(* Out[30]= SeriesData[x, 0, {1, c[3], c[4], c[5], c[6], c[7]}, 2, 8, 1] *)

What InverseSeries does, as I recall, is some form of Newton's method, effectively doubling the number of terms at each step. My guess as to bottleneck is it is probably trying for internal "simplifications" e.g. using Together, that may work often but do not help much in the case of independent symbolic coefficients. But this is just a guess, I have not tested it out.

As for having a "general form" applicable to the case of independent symbolic coefficients, it would not surprise me if a direct method applies, and is substantially more efficient. The "magic" method could be such an approach. MathWorld has an entry containing useful references for this.

POSTED BY: Daniel Lichtblau

Now that I am just talking about numbers and not drawing any pictures, I think this thread may be losing interest. Also we haven't heard from the relevant experts as to why / how I am out-competing a usually-strong Mma function. But I talked breifly with Peter Luschny who found an interesting connection to the Pochammer function. Working from there and searching the OEIS more I obtained the following closed form in terms of MultiFactorial:

MultiFactorial[n_, nDim_] := 
Times[n, If[n - nDim > 1, MultiFactorial[n - nDim, nDim], 1] ]

AlternateGeneralT[n_, nLead_] := 
Table[MultiFactorial[i + nLead j + 2, nLead]/
MultiFactorial[i + 2, nLead], {i, 0, n}, {j, 0, i}] // TableForm

And the following Prints:

Grid[Partition[(Abs@GeneralT[5, #] // TableForm) & /@ Range[12], 4], 
 Frame -> All]

Grid[Partition[AlternateGeneralT[4, #] & /@ Range[12], 4], 
 Frame -> All]

Yield the same triangles

Number Triangles

Thank you Wolfram Community Moderation Team, for showing me a little bit of respect.

Conclusion

To reiterate my points:

  • What is the Mathematica Algorithm?
  • Why is Mathematica going so slow?
  • Can we change the Mma Source Code to use one of my definitions or something similar?
  • Why aren't any of the triangles above already in OEIS?
  • Is this reduction known?

I find the discovery enjoyable, and look forward to finding out just why it works this way. More work ahead . . .

Next, I'm going back to physics research, where I apply these triangles of numbers to perturbation calculations of planetary motion. Furthermore, I have another idea for applying series inversion to make it easy to compute integrals along Rotational Energy Surfaces.

POSTED BY: Brad Klee

I can comment on the speed issue, or, more accurately, on the lack of prior comment about the speed issue. This is the sort of thing we are able to investigate in principle, but it absolutely requires clear concise examples. What I see is code from a few posts, some of it seemingly special-case for certain series (I'm not sure, but that's what it looks like), some definitions that seem incomplete. It's quite difficult to dive into this, even if the code produces nice graphics and shows good speed. In particular, a small example that compares two methods, where neither is special-cased for a particular series, would be a minimal requirement.

I will say that it is entirely possible that InverseSeries suffers from a deficiency in efficiency, so to speak. But again, any investigation will require a short concise example, where the desired outcome is made clear. This can be done perhaps by showing a small case for which InverseSeries performance is fine, and simply indicating what size parameter needs to be adjusted in order for the speed issue to manifest.

POSTED BY: Daniel Lichtblau

HAKMEM: General Algorithms for Leading Term $x^n$

Series

$$ y = x^n + c_{n+1}\; x^{n+1}+ c_{n+2} \;x^{n+2} + c_{n+3}\; x^{n+3} + \ldots $$

Mathematica Built In

MmaExpand[m_, nLead_] := Expand[Refine[Normal[
     InverseSeries[ Series[f[x], {x, 0, m + nLead}] /. Join[
        {f[0] -> 0, (D[f[x], {x, nLead}] /. x -> 0) -> (nLead!)},
        (D[f[x], {x, #}] /. x -> 0) -> 0 & /@ Range[1, nLead - 1],
        ((D[f[x], {x, #}] /. x -> 0) -> #! c[#]) & /@ 
         Range[nLead + 1, m + nLead]]]] /. x -> y^nLead, y > 0]]

Conceptually Simple Multinomial Method

RExp[n_] := Expand[y Plus[R[0], Total[y^# R[#] & /@ Range[n]]]]

RCalc[0, _] = 1;
RCalc[n_, nLead_] := With[{basis = Subtract[Tally[Join[Range[n + nLead], #]][[All, 2]],
Table[1, {n + nLead}]] & /@ Select[IntegerPartitions[n + nLead], Length[#] > nLead - 1 &][[2 ;; -1]]},
Total@ReplaceAll[Times[-1/nLead, Multinomial @@ #, c[Total[#]],
Times @@ Power[RSet[# - 1] & /@ Range[n + nLead], #]] & /@     basis, {c[nLead] -> 1}]]

MultinomialExpand[n_, nLead_] := ( Clear@RSet; Set[RSet[#], Expand@RCalc[#, nLead]] & /@ Range[0, n];
Expand[RExp[n] /. R[m_] :> RCalc[m, nLead]])

Frobenius-Stirling Conjectural Magic

Also see Peter Luschny, Stirling-Frobenius Cycle Numbers.

SF[0, 0, m_] = 1;
SF[n_, k_, m_] := SF[n, k, m] = If[Or[k > n, k < 0], 0,
If[And[n == 0, k == 0], 1,
SF[n - 1, k - 1, m] + (m*n - 1)*SF[n - 1, k, m]
]]

GeneralPolys[n_, nLead_] := MapThread[Dot, {
Table[SF[i, j, nLead], {i, 0, n - 1}, {j, 0, i}],
Table[ x^j, {i, 0, n - 1}, {j, 0, i}]}]

GeneralT[n_, nLead_] := With[{ps = GeneralPolys[n, nLead]},
Table[(-1)^j  ps[[j]] /. {x -> i + 2}, {i, 1, n}, {j, 1, i} ]
]

Basis[n_, nLead_] := Total /@ Map[
Times[nLead^(-Length[#]) 1/Times @@ (Tally[#][[All, 2]]!),
Times @@ (# /. x_Integer :> c[x + nLead])] &, Function[{a},
Select[IntegerPartitions[n], Length[#] == a &]] /@ Range[n], {2}]

MagicExpand[n_, nLead_] := Expand[y + Dot[MapThread[Dot,
{GeneralT[n, nLead], Basis[#, nLead] & /@ Range[n]}], 
y^Range[2, n + 1]]];

Compare

In[]:= ColumnForm[MmaExpand[2, #] /. y -> (# y) & /@ Range[5]]
ColumnForm[MultinomialExpand[2, #] /. y -> (# y) & /@ Range[5]]
ColumnForm[MagicExpand[2, #] /. y -> (# y) & /@ Range[5]]

Out[]= ColumnForm[{
y - y^2 c[2] + 2 y^3 c[2]^2 - y^3 c[3], 
2 y - 2 y^2 c[3] + 5 y^3 c[3]^2 - 4 y^3 c[4], 
3 y - 3 y^2 c[4] + 9 y^3 c[4]^2 - 9 y^3 c[5], 
4 y - 4 y^2 c[5] + 14 y^3 c[5]^2 - 16 y^3 c[6], 
5 y - 5 y^2 c[6] + 20 y^3 c[6]^2 - 25 y^3 c[7]}]

Out[]= ColumnForm[{
y - y^2 c[2] + 2 y^3 c[2]^2 - y^3 c[3], 
2 y - 2 y^2 c[3] + 5 y^3 c[3]^2 - 4 y^3 c[4], 
3 y - 3 y^2 c[4] + 9 y^3 c[4]^2 - 9 y^3 c[5], 
4 y - 4 y^2 c[5] + 14 y^3 c[5]^2 - 16 y^3 c[6], 
5 y - 5 y^2 c[6] + 20 y^3 c[6]^2 - 25 y^3 c[7]}]

Out[]= ColumnForm[{
y - y^2 c[2] + 2 y^3 c[2]^2 - y^3 c[3], 
2 y - 2 y^2 c[3] + 5 y^3 c[3]^2 - 4 y^3 c[4], 
3 y - 3 y^2 c[4] + 9 y^3 c[4]^2 - 9 y^3 c[5], 
4 y - 4 y^2 c[5] + 14 y^3 c[5]^2 - 16 y^3 c[6], 
5 y - 5 y^2 c[6] + 20 y^3 c[6]^2 - 25 y^3 c[7]}]

In[]:= AbsoluteTiming[SameQ[
MagicExpand[15, #] & /@ Range[5],
MultinomialExpand[15, #] & /@ Range[5]
]]

Out[]= {4.61176, True}

In[]:= AbsoluteTiming[SameQ[
MmaExpand[10, #] & /@ Range[5],
MultinomialExpand[10, #] & /@ Range[5]
]]

Out[]= {9.79921, True}

Timing test for $n=1,2,3,4,5$

Time test

Time Vs. order of magnitude.

Not only is the mystery method the fastest, apparently it is also the most stable under change of leading term. I may try to prove this Stirling-Frobenius connection in my dissertation, but before I get into that I'm wondering: Is this algorithm already proven somewhere? Did Ferdinand Frobenius know about this? Are there Mma experts around here with an answer? Maybe Daniel Lichtblau will care to comment? Considering that Mma is going relatively slow and the optimization code is easy, I'm thinking that there's a possibility the methods here are not well known, possibly new.

POSTED BY: Brad Klee

Adapted Multinomial Algorithm

Comparing $n=2$ case above, and $n=1$ case here, it seems like the following simple algorithm could easily be generalized to revert any series with leading term order $n$ :

RExp[n_] := Expand[b Plus[R[0], Total[b^# R[#] & /@ Range[n]]]]
RCalc[n_] := 
With[{basis = Subtract[Tally[Join[Range[n + 1], #]][[All, 2]],
Table[1, {n + 1}]] & /@ IntegerPartitions[n + 1][[2 ;; -1]]},
Total@(Times[-1, Multinomial @@ #, v[Total[#]],
Times @@ Power[RSet[# - 1] & /@ Range[n + 1], #]] & /@ basis)]
InvSeries[n_] := AbsoluteTiming[
Clear@RSet; RSet[0] = 1; 
Set[RSet[#], Expand@RCalc[#]] & /@ Range[n];
Expand[RExp[n] /. R -> RSet]]

Updated Time Test

Time Test

Conceptually Simple method using multinomial coefficients fits in between slower Mma code and faster unproven method. This is the same time-behavior as seen in the case with leading quadratic term.

Challenge: Calculate series reversion with leading cubic leading term. Is there a regularization of the partition structure? Assuming yes, what is it? How does it relate to results for $n=1,2$ ?

POSTED BY: Brad Klee

Lifting Megabytes in the Gymnasium Pt. 2

Confer: A283298, A283247.

Stronger Defintion for Symbolic Reversion of Power Series ( conjecture )

By analogy to the above, we rework the case with leading linear term, as presented on Mathworld, and in Abramowitz and Stegun.

PolysA257635[n_] :=  Expand[-#! Binomial[2 # - x, #] & /@ Range[0, n - 1]]
TA283298[n_] := With[{polys = PolysA257635[n]},   Table[polys[[j]] /. x -> 2 j + i, {i, 1, n}, {j, 1, i}]]

Functions for Comparison

InvExp = InverseSeries[ Series[f[x], {x, 0, 11}] /. (D[f[x], {x, n_}] /. x -> 0) :> (n! \[Epsilon][n]) /.
{  f[0] -> 0, \[Epsilon][1] -> 1}];

RingGens[n_] := Times @@ (\[Epsilon] /@ #) & /@ (IntegerPartitions[n] /. x_Integer :> x + 1)

Basis[n_] := 
Total /@ Map[
Times[(*2^(Length[#]-1)*)1/Times @@ (Tally[#][[All, 2]]!), 
Times @@ (# /. x_Integer :> \[Epsilon][x + 1])] &, 
Function[{a}, Select[IntegerPartitions[n], Length[#] == a &]] /@ 
Range[n], {2}]

A111785[n_] := Expand[x + Dot[
MapThread[Dot, {TA283298[n], Basis[#] & /@ Range[n]}],
x^Range[2, n + 1]     ]];

A111785T[n_] :=  With[{exp = A111785[n]},  
Function[{a},    Coefficient[Coefficient[exp, x, a + 1], #] & /@ RingGens[a]] /@  Range[n]]

Brute Force

In[] = SameQ[A111785[10], Expand[Normal[InvExp]]]
Out[] = True

Application to Arctangent

\[Epsilon]List[
n_] := \[Epsilon][#] -> 
Coefficient[Normal@Series[Tan[x], {x, 0, n + 1}], x, #] & /@ 
Range[n + 1];

fList = A111785[#] /. \[Epsilon]List[#] & /@ (2 Range[15]);

Colors = Blend[{Yellow, Orange}, #/15] & /@ Range[15];

Show[
Plot[ArcTan[x], {x, -Pi/2, Pi/2}, PlotStyle -> Red],
MapThread[
Plot[#1, {x, -Pi/2, Pi/2}, PlotStyle -> #2] &, {fList, Colors}],
ImageSize -> 800
]

Arctangent

Time Test and ByteCount

Time Test

The following command practically hangs on my personal computer:

In[] = AbsoluteTiming[ N[ByteCount[    Expand[Normal[
InverseSeries[       Series[f[x], {x, 0, 15}] /. (D[f[x], {x, n_}]
/. x -> 0) :> (n! \[Epsilon][n]) /. {     f[0] -> 0, \[Epsilon][1] -> 1}]]]]/10^6]]
Out[]= $Aborted

How about the strong man function, how many Megabytes?

In[]= AbsoluteTiming[ N[ByteCount[A111785[35]]/10^6]]
Out[]= {6.73205, 47.5168}

Almost 50 in under 10s. Now that's a Powerful Power series! Conjectural magic and mysterious binomial number theory wins again !

POSTED BY: Brad Klee

enter image description here - Congratulations! This post is now Staff Pick! Thank you for your wonderful contributions. Please, keep them coming!

POSTED BY: Moderation Team

Timing Test

Mysterious Mathematica Function ( ? )

Fast[m_] := Expand[Refine[Normal[
InverseSeries[ Series[f[x], {x, 0, m + 2}] /. Join[
{f[0] -> 0, (D[f[x], x] /. x -> 0) -> 
0, (D[f[x], {x, 2}] /. x -> 0) -> 
1/2}, ((D[f[x], {x, #}] /. x -> 0) -> #! \[Epsilon][#]) & /@
Range[3, m + 2]]]] /. x -> b^2/4, b > 0]]

Multinomial Method, A276738

(* with functions above *)
Faster[n_] := Module[{},
Clear@RSet; RSet[0] = 1; Set[RSet[#], Expand@RCalc[#]] & /@ Range[n];
Expand[b + RExp[n] /. R -> RCalc /. v -> \[Epsilon]]
]

Mystery via Triangle A028338 ( conjecture )

A028338[n_, k_] := Sum[(-2)^(n - i) Binomial[i, k] StirlingS1[n, i], {i, k, n}]

Polys[n_] := Total /@ Table[
A028338[i, j] (-1)^(j + 1) x^(j), {i, 0, n - 1}, {j, 0, i}];

T2[n_] := With[ {polys = Polys[n]},
Table[polys[[j]] /. x -> (2 j + i), {i, 1, n}, {j, 1, i}] ]

Basis[n_] := Total /@ Map[
Times[2^(Length[#] - 1)/Times @@ (Tally[#][[All, 2]]!), 
Times @@ (# /. x_Integer :> \[Epsilon][x + 2])] &,
Function[{a}, Select[IntegerPartitions[n], Length[#] == a &]] /@ 
Range[n],{2}]

Fastest[n_] := Expand[Plus[b,
Total@With[{tab = T2[n]}, 
2 Expand[Basis[#].tab[[#]] b^(# + 1)] & /@ Range[n]]]]

Output

In[14]:= Fast[2]
Faster[2]
Fastest[2]

Out[14]= b - 2 b^2 \[Epsilon][3] + 10 b^3 \[Epsilon][3]^2 - 
2 b^3 \[Epsilon][4]

Out[15]= b - 2 b^2 \[Epsilon][3] + 10 b^3 \[Epsilon][3]^2 - 
2 b^3 \[Epsilon][4]

Out[16]= b - 2 b^2 \[Epsilon][3] + 10 b^3 \[Epsilon][3]^2 - 
2 b^3 \[Epsilon][4]

More Testing

In[17]:= SameQ[Fast[5], Faster[5]]
SameQ[Faster[10], Fastest[10]]

Out[17]= True

Out[18]= True

Racing Functions

tFast = AbsoluteTiming[
Fast[#]
][[1]] & /@ Range[11]

tFaster = AbsoluteTiming[
Faster[#]
][[1]] & /@ Range[22]

tFastest = AbsoluteTiming[
Fastest[#]
][[1]] & /@ Range[35];

ListLinePlot[{tFast, tFaster, tFastest}, PlotRange -> All, 
ImageSize -> 600]

Race

So we see that $ T_{Mma} \gg T_{A276738} \gg T_{A028338} \approx 2(sec) $, where $T$ is the time to expand $30$ orders of magnitude !

Looks like it could be possible to improve the timing of "InverseSeries", at least in this interesting and useful case. It's also worthwhile to note that my code here for the mystery method could still be slow as it is easy to generate coefficients for $100$ orders of magnitude:

In[65]:= AbsoluteTiming[T2[100]][[1]]
Out[65]= 0.896722  

Wow! That is really fast.

I've seen GOSPER around here talking about fractals again. I'm wondering if we will hear from him regarding A028338, which he originally authored. For what purpose, Bill ?

POSTED BY: Brad Klee

Deformation Animation

Deformation Animation

Time animation condition: $$ \frac{d}{dt} \Psi \approx \pm c_1, $$ along direction of maximum deformation, $(1,1,1)$. This requires, $$ \frac{d}{dt} E \neq c_2, $$
with arbitrary constants $ c_1, c_2 $.

Color Gradient Surface Tilings

Colored Tiling

Inspired by Old Veteran Paul Klee Einsame Blute .

POSTED BY: Brad Klee
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract