# [GIF] Plotting the Contours of Deformed Hyperspheres

GROUPS:

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}}]},
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];

## Octahedral

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

1 year ago
13 Replies

## 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$.

Inspired by Old Veteran Paul Klee Einsame Blute .

1 year ago

# 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]

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 ?

1 year ago
 - Congratulations! This post is now Staff Pick! Thank you for your wonderful contributions. Please, keep them coming!
1 year ago

# 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],
Plot[#1, {x, -Pi/2, Pi/2}, PlotStyle -> #2] &, {fList, Colors}],
ImageSize -> 800
]

## Time Test and ByteCount

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]]