Have I waited long enough? Because I'm still hoping to find the C∞ solution, and I did some more work on this myself.dummy_index wrote: ↑Sat Aug 27, 2022 10:08 pmPlease wait a while as I would like to derive the C3 class solution on my own.
I figured out how to code Faà di Bruno's formula in Maxima:
Code: Select all
halfangles: true$
g[n] := at(diff(asin(cos(a)/sqrt(2)), a, n), a=acos(1/3)/2)$
bell[n,k] := if k=1 then g[n] else fullratsimp(sum(binomial(n,j) * g[j] * bell[n-j,k-1], j, 1, n-1) / k)$
f[n] := fullratsimp(sum(f[k] * bell[n,k], k, 1, n-1) / (1 - bell[n,n]))$ /* Note that this solves f(a) = f(g(a)), given g(a) = a, ignoring constant and linear terms in the actual formula we want to solve, and so is valid only for n >= 2. n = 1 is manually patched in below. */
f[1]: %pi/acos(1/3) / 3$
f[1] = %pi/acos(1/3) * 1/3
f[2] = %pi/acos(1/3) * -1/3 / sqrt(2)
f[3] = %pi/acos(1/3) * 1/9
f[4] = %pi/acos(1/3) * -1/5 / sqrt(2)
f[5] = %pi/acos(1/3) * -29/9
f[6] = %pi/acos(1/3) * 439/77 / sqrt(2)
f[7] = %pi/acos(1/3) * -42979/4257
f[8] = %pi/acos(1/3) * -65951/3655 / sqrt(2)
f[9] = %pi/acos(1/3) * 883085651/1375011
Unfortunately these follow no obvious pattern, although they do appear to stay rational multiples of either %pi/acos(1/3) (odd derivatives) or %pi/acos(1/3)/sqrt(2) (even derivatives).
Furthermore this computes the derivatives at a = acos(1/3)/2, whereas a Taylor series would be much simpler (containing only odd-exponent terms) around a = %pi/4. Plus, given the non-analytic nature of the function, the Taylor series around a = acos(1/3)/2 probably doesn't actually converge.
Suppose we try to define f (more accurately, f_2 from before, valid between acos(1/3)/2 and pi/2-acos(1/3)/2) as
f(a) := sum(c[m] / (2*m-1)! * (a-%pi/4)^(2*m-1), m, 1, inf) + %pi/4
Then c[m] need to satisfy the following relationship relative to the previously-computed f[n]:
f[n] = sum(c[m] / (2*m-1-n)! * (acos(1/3)/2-%pi/4)^(2*m-1-n), m, ceiling((n+1)/2), inf)
Actually solving this to infinity is difficult, but if we limit the summation to, say, five terms, then we get:
c[1] = -(%pi*(145*2^(9/2)*acos(1/3)^4-145*2^(11/2)*%pi*acos(1/3)^3-44352*acos(1/3)^3+435*2^(7/2)*%pi^2*acos(1/3)^2+66528*%pi*acos(1/3)^2-4785*2^(13/2)*acos(1/3)^2-145*2^(7/2)*%pi^3*acos(1/3)-33264*%pi^2*acos(1/3)+4785*2^(13/2)*%pi*acos(1/3)-5892480*acos(1/3)+145*sqrt(2)*%pi^4+5544*%pi^3-4785*2^(9/2)*%pi^2+2946240*%pi-495*2^(31/2)))/(1485*2^(31/2)*acos(1/3))
c[2] = (%pi*(1160*acos(1/3)^3-1740*%pi*acos(1/3)^2-297*2^(11/2)*acos(1/3)^2+870*%pi^2*acos(1/3)+297*2^(11/2)*%pi*acos(1/3)-100320*acos(1/3)-145*%pi^3-297*2^(7/2)*%pi^2+50160*%pi-17325*2^(11/2)))/(380160*acos(1/3)*(2*acos(1/3)-%pi))
c[3] = -(%pi*(232*acos(1/3)^3-348*%pi*acos(1/3)^2-99*2^(9/2)*acos(1/3)^2+174*%pi^2*acos(1/3)+99*2^(9/2)*%pi*acos(1/3)-12320*acos(1/3)-29*%pi^3-99*2^(5/2)*%pi^2+6160*%pi-1155*2^(11/2)))/(264*acos(1/3)*(2*acos(1/3)-%pi)^3)
c[4] = (8*%pi*(1160*acos(1/3)^3-1740*%pi*acos(1/3)^2-99*2^(13/2)*acos(1/3)^2+870*%pi^2*acos(1/3)+99*2^(13/2)*%pi*acos(1/3)-36960*acos(1/3)-145*%pi^3-99*2^(9/2)*%pi^2+18480*%pi-3465*2^(11/2)))/(33*acos(1/3)*(2*acos(1/3)-%pi)^5)
c[5] = -(1792*%pi*(1160*acos(1/3)^3-1740*%pi*acos(1/3)^2-297*2^(9/2)*acos(1/3)^2+870*%pi^2*acos(1/3)+297*2^(9/2)*%pi*acos(1/3)-26400*acos(1/3)-145*%pi^3-297*2^(5/2)*%pi^2+13200*%pi-2475*2^(11/2)))/(33*acos(1/3)*(2*acos(1/3)-%pi)^7)
...Yeah. Numerically, this gives:
c[1] = 0.778240456947945
c[2] = 7.423254731907882
c[3] = -1488.676285956136
c[4] = 619910.1478436295
c[5] = -214811061.0840599
Compare the earlier C2 version, which had:
c[1] = 0.799609546057049
c[2] = 3.540212531753472
Or the even earlier C1 version, which had:
c[1] = 0.8507165520199232
Note that c[m] grows even faster than (2*m-1)! does, so even if we redefined the Taylor summation to omit that part, they'd still get pretty high.
Going the other way, a C20 version has:
c[1] = 0.764495900777279 = 0.764495900777279 * 1!
c[2] = 16.2937621351 = 2.715627022516667 = 3!
c[3] = -20022.43327616207 = -166.8536106346839 * 5!
c[4] = 7.076072603208672E+7 = 14039.82659366799 * 7!
c[5] = -3.92158350587567E+11 = -1080683.285349336 * 9!
c[6] = 2.852412677766556E+15 = 71458951.56341581 * 11!
c[7] = -2.489653753710991E+19 = -3998145876.935229 * 13!
c[8] = 2.466375940530462E+23 = 188607806414.5754 * 15!
c[9] = -2.650236337102935E+27 = -7451026175678.147 * 17!
c[10] = 3.14376953136687E+31 = 258437826168181.4 * 19!
c[11] = -4.142165954208365E+35 = -8107437009650615 * 21!
c[12] = 8.454279234842538E+36 = 327025907504003 * 23!
c[13] = -1.603867150320828E+39 = -103400517808757.8 * 25!
c[14] = 1.686363403437531E+41 = 15487038494825.02 * 27!
c[15] = -6.217355432072262E+44 = -70318059188591.4 * 29!
c[16] = -1.637231301705888E+45 = -199107798481.9793 * 31!
c[17] = 1.401911494166907E+48 = 161448832774.4398 * 33!
c[18] = 9.29084015974988E+49 = 8991296930.976983 * 35!
c[19] = -3.282134015861778E+51 = -238462139.948875 * 37!
c[20] = -8.760449740820171E+52 = -4294783.99078278 * 39!
This does feel like it converges to something, but...
Oh, the command for computing those:
Code: Select all
coefficients(continuity) := solve(makelist(f[n] = sum(c[m] / (2*m-1-n)! * (acos(1/3)/2-%pi/4)^(2*m-1-n), m, ceiling((n+1)/2), continuity), n, 1, continuity), makelist(c[n], n, 1, continuity))[1]$