RE: Re: RE: RE: Functional decomposition (solving f[f[x]] = g[x] for given g)
- To: mathgroup at smc.vnet.net
- Subject: [mg71972] RE: [mg71885] Re: [mg71872] RE: [mg71839] RE: [mg71764] Functional decomposition (solving f[f[x]] = g[x] for given g)
- From: "Ingolf Dahl" <ingolf.dahl at telia.com>
- Date: Wed, 6 Dec 2006 06:04:05 -0500 (EST)
- Organization: Göteborg University
- Reply-to: <ingolf.dahl at telia.com>
Hi again dear MathGroup,
I have checked the suggested procedure below for solving the equation
f(f(x))=exp(x), and it works smoothly.
Maybe this is already well-known. Maybe it is known, but not implemented in
Mathematica code. Maybe it is included in the book by Aczel, I have not been
able to check that yet, since I have not the book available. If you know
where to find an already published solution, please inform me.
But anyway, I have not spent too many hours on it, and it is a nice
demonstration of the powers of Mathematica. Thus it is relevant to publish
it here.
Moreover, the solution I obtain is characterized by some mathematical
constants, which should be thumbprints for uniqueness. I have not been able
to localize these constants anywhere else.
But what I have found is that the solution is related to half-integer values
of "tetration" in Wikipedia and PowerTower in MathWorld (naturally extending
into other fractional values). Wikipedia states; "At this time there is no
commonly accepted solution to the general problem of extending tetration to
the real or complex numbers, although it is an active area of research."
What I claim is that (in the Mathematica notation from FurtherExamples to
"Power" in Help)
PowerTower[E,-1/2] = "HalfwayFromZero" =
0.4987433645316708791828373750416860523968424`43.
PowerTower[E,-3/2] = "HalfwayToZero" =
-0.6956636150637165750917227913152630811483903`42.
(PowerTower[E,-2] = -Infinity, PowerTower[E,-1] = 0, PowerTower[E,0] = 1 by
natural extension)
Here is my solution:
I started by doing a routine for decomposition of some functions g(x) into q
equal functions. Thus f(f(f...(f(x))...))=g(x), with q f:s. First the series
expansion to the power n is derived. Then the function is composed p times
with itself. We thus get the series expansion of the function to the
decompositional power p/q. The requirement on g is that g[0)=0, with
positive derivative at origin and with a Maclaurin expansion around origin.
Attributes[FractionalSeries]= {HoldFirst};
FractionalSeries[ func_/; And[ func[0.]==0., SeriesCoefficient[ Series[
func[var], { var,0,1}],1]>0],x_,
{ p_Integer/; p>0, q_Integer/; q>0}, n_Integer/; n>0]:=
Module[ { a,f1,fn,f1s},
f1[var_]= Sum[ var^i* a[i], { i,1,n}];
fn[var_]:=
Nest[ f1, var,q];
a[1]= a[1] /. Solve[ Reduce[ { SeriesCoefficient[ Series[ fn[var], {
var,0,1}],1]== SeriesCoefficient[ Series[ func[var], { var,0,1}],1],
a[1]>0}, a[1]], a[1]][ [1]];
Do[ a[i]= a[i] /. Solve[ { SeriesCoefficient[ Series[ fn[var], {
var,0,i}],i]== SeriesCoefficient[ Series[ func[var], { var,0,i}],i]}, a[i]][
[1]], { i,2,n}];
Series[ f1[x], { x,0,n}];
ComposeSeries@@ Table[ Series[ f1[x], { x,0,n}], {p}]
Example of use: Sinh is such a function, and let us take it to the
decompositional power 2/3. Define
ftwothirds[y_] := Normal[FractionalSeries[Sinh[#1] & , var, {2, 3}, 29]] /.
{var -> y};
Then
ftwothirds[ftwothirds[ftwothirds[0.3]]] - Sinh[Sinh[0.3]]
evaluates to
5.551115123125783*^-17
and
ftwothirds[ftwothirds[ftwothirds[SetAccuracy[0.3, 60]]]] -
Sinh[Sinh[SetAccuracy[0.3, 60]]]
evaluates approx to
1.65*^-24
The problem is that Exp is not such a function, since Exp[0} == 1. But Exp
behaves asymptotically as 2*Sinh. Let us find the series expansion of 2*Sinh
to the decompositional power 1/2. (There is a more complicated relation
between the series expansions for 2*Sinh and Sinh than I thought from the
beginning.)
Simplify[FractionalSeries[2*Sinh[#1] & , y, {1, 2}, 29]]
(returns long expression)
If we use this expression, we might build a function f2sinh to extend the
series expansion to whole real axis. f2sinh should now be THE smooth
solution to the equation f2sinh[f2sinh[x]] == Sinh[x] .
f2sinh[x_Real] := Module[{y, steps},
y = If[MachineNumberQ[x], SetAccuracy[x, 24], x]; steps = 0;
While[Abs[y] > 0.07, steps = steps + 1; y = ArcSinh[y/2]];
y = fpoly[y];
(If[MachineNumberQ[x], N[#1], #1] & )[Nest[2*Sinh[#1] & , y, steps]]]
where
fpoly[y_] := Sqrt[2]*y + y^3/(9*Sqrt[2]) - (11*y^5)/(2700*Sqrt[2]) +
(479*y^7)/(1020600*Sqrt[2]) - (426679*y^9)/(6246072000*Sqrt[2]) +
(85589473*y^11)/(7557747120000*Sqrt[2]) -
(1403953837837*y^13)/(689720002171200000*Sqrt[2]) +
(1031743413846247*y^15)/(2669216408402544000000*Sqrt[2]) -
(1496209724379058102331*y^17)/(19591834900362000756480000000*Sqrt[2]) +
(17812978070670318115902941*y^19)/(1145769688642970528240463360000000*Sqrt[2
]) -
(7533260191446785723514332836231*y^21)/(232533958310090868706402038912000000
0000*Sqrt[2]) +
(136561976266545063524556174223350449*y^23)/(1985445516947137017537762080849
52576000000000*Sqrt[2]) -
(719728529960035731908565151209614126425839277*y^25)/(4853778747825775229320
514450035010550004224000000000000*Sqrt[2]) +
(100427683629367623630378641852188624798406912708777*y^27)/(3101826723913052
963398192041352673632022399364096000000000000*Sqrt[2]) -
(230559003784630013844879523283194881704293155047113459273*y^29)/(3226456204
1869955579438812968736894041769959385166794752000000000000*Sqrt[2]) +
SetAccuracy[0., 44]*(y/SetAccuracy[0.07, 60])^29;
The limit 0.07 in this definition is a bit arbitrarily chosen, and you might
decrease it to increase accuracy, which anyway is good along the whole real
axis. I have estimated the truncation error to be as large as the last term
in the series.
Now we can assume that Exp to a fractional decompositional power
asymptotically should be identical to 2*Sinh to the same fractional
decompositional power, with a relative error that is less than the relative
error of Exp compared to 2*Sinh, if the decompositional power is less than
one. Thus for any x, we might climb up to a large value of Exp...[Exp[x]]
(say x>60), fetch a value using f2sinh[x], and then climb down again. Thus
we might calculate the (smooth) solution to the equation fexp[fexp[x]]==
Exp[x] as
fexp[x_Real] :=
Module[{y, steps},
y = If[MachineNumberQ[x], SetAccuracy[x, 24], x]; steps = 0;
While[y < 60, steps = steps + 1; y = Exp[y]];
y = f2sinh[y]*(1 + SetAccuracy[0., y/1.15]);
(If[MachineNumberQ[x], N[#1], #1] & )[Nest[Log[#1] & , y, steps]]]
The limit of fexp[x] when x->-Infinity is the constant HalfwayToZero defined
above.
And fexp[HalfwayToZero] = 0,
fexp[0]= HalfwayFromZero,
fexp[HalfwayFromZero] = 1,
fexp[1] = Exp[HalfwayFromZero],
fexp[Exp[HalfwayFromZero] = E und so weiter.
Any devils hiding anywhere?
One little devil is that fexp[x] approaches HalfwayToZero fast when z ->
-Infinity, so there easily is a loss of accuracy. For large negative values
of x, fexp will approximately follow the formula
fexp[x] ~ HalfwayToZero + Exp[x]/f´[HalfwayToZero]
where
f´[HalfwayToZero] = 0.5688037784984972069377219276490924736217`40.
For calculations of derivatives etc. of f(x) in the region around x = 0, a
slightly simpler formula might be used instead of f(x):
fexpsimpl[x_Real] := Nest[Log, Nest[2Sinh[#] &,
fpoly[Nest[ArcSinh[#/2] &, Nest[Exp, x, 4], 8]], 8], 4]
Thus we can obtain the zeroth, first second and third derivatives of fexp(x)
at x=0:
{0.498743364531670879182837375041686052396842401697695,
0.87682850111901739234161348908444926737095344685410,
0.489578037438247862041903165991088781934758253721,
0.1222616742183662064640316453590024071635142693}
It should be possible to calculate Exp to other fractional decompositional
powers, and to calculate PowerTower for other first arguments than E and for
fractional second argument, following the general scheme here.
That all for today!
Ingolf Dahl
> -----Original Message-----
> From: Gregory Duncan [mailto:gmduncan at gmail.com]
> Sent: den 2 december 2006 11:11
> To: mathgroup at smc.vnet.net
> Subject: [mg71885] Re: [mg71872] RE: [mg71839] RE: [mg71764]
> Functional decomposition (solving f[f[x]] = g[x] for given g)
>
> The problem falls generally into the the area of functional equations.
> See J. Aczel
> Functional Equations: History, Applications and Theory or any
> of his other books. One can also try to find a fixed point in
> a function space.
>
>
> On 12/1/06, Ingolf Dahl <ingolf.dahl at telia.com> wrote:
> > Some further comments on the equation f(f(x)) = g(x):
> > For the case g(x) = exp(x) some work have been done, see
> > http://www.math.niu.edu/~rusin/known-math/index/39-XX.html ,
> > http://www.math.niu.edu/~rusin/known-math/99/sqrt_exp
> > and
> http://www.newton.dep.anl.gov/newton/askasci/1993/math/MATH023.HTM .
> > (found by a Google search for "f(f(x))")
> >
> > However, I have not penetrated deeply enough into all references to
> > really see if there are good methods to find (precise) numerical
> > values for at least one solution of the equation. Or a plot...
> >
> > But some thoughts further, somewhat fishing on deep water:
> If we have
> > a function g(x) that is such that g(0)=0 (or g(x0)=x0 for
> some x0) and
> > where
> > g(x) also has a Maclaurin expansion with a positive first
> coefficient,
> > it should be easy to find a Maclaurin series which composed
> to itself
> > is identical to that of g(x). Examples of such functions g(x) are
> > sin(x) and sinh(x). Then there is the question if this series
> > converges in some range, I have not tested yet. (In the links above
> > they refer to a theorem by Polya, stating that convergence
> range must
> > be limited.) For some other functions, we could also think
> of expansions in x^(Sqrt(2)).
> > Suppose that we find such a function f(x) for g(x) =
> sinh(x), name it
> > fsinh(x), with convergence in some region near origo. Then
> it might be
> > possible to extend that to the whole axis if we for any x
> apply ginv
> > (the inverse of g) a number of times until we reach the convergence
> > region. (we might have to use extended precision.) There we
> apply f,
> > and then apply g the same number of times as we did with
> ginv. Then we
> > should have obtained f(x).
> >
> > If we now look at the case g(x) = exp(x), we can note that
> if x > 20,
> > exp and 2*sinh are essentially identical (to machine
> precision). Thus
> > we could assume that fsinh(x)*Sqrt[2] is a good
> approximation of the
> > solution f(x) for g(x)=exp(x) and x>20, and then we might
> use the same
> > procedure as for
> > g(x)=x^2+1 to extend this to zero. With the power of Mathematica
> > available, just a few lines of code maybe are all that is
> needed. But
> > maybe some devils are hiding in the details?
> >
> > Is there anyone who knows more about this?
> >
> > Ingolf Dahl
> >
> >
> > > -----Original Message-----
> > > From: Ingolf Dahl [mailto:ingolf.dahl at telia.com]
> > > Sent: den 30 november 2006 12:06
> > > To: mathgroup at smc.vnet.net
> > > Subject: [mg71839] RE: [mg71764] Functional decomposition
> (solving
> > > f[f[x]] = g[x] for given g)
> > >
> > > Hi Kelly,
> > > Here is a rapid amateur answer (but without scary rhymes,
> as in my
> > > last answer). If you are satisfied with the numerical
> values of the
> > > function, you can investigate
> > >
> > > g[x_]:=1+x^2;
> > > ginv[x_]:=Sqrt[x-1];
> > > f[x_Real]:=Nest[ginv,Power[Nest[g,x,9],Sqrt[2]],9]
> > >
> > > f[x_Real] appears to be a good solution to f[f[x]] == g[x] to
> > > approximate machine precision. The solution is based on two
> > > observations:
> > > 1) Nest[g,x,9] is equivalent to Nest[f,x,18], and we reach very
> > > large values when this is calculated. There g[x] is approximately
> > > equal to x^2 (to machine precision), and f[x] can be
> approximated by
> > > x^Sqrt[2]. Thus we might reach a good approximation of
> Nest[f,x,19]
> > > by Power[Nest[g,x,9],Sqrt[2]].
> > > 2) We can then apply the inverse function ginv of g 9
> times to reach
> > > f[x].
> > > The relative error is reduced by the application of the inverse
> > > function.
> > >
> > > We can note that f[0.] is 0.6420945043908285` which do
> not appear to
> > > be a known mathematical constant. But now it is.
> > > (If someone can check... I cannot reach the "The Inverse Symbolic
> > > Calculator" at CESM. Is it down?)
> > >
> > > Generalize and write f[x] as DecompositionPower[x,1/2] and the
> > > solution to f[f[f[x]]]\[Equal]g[x] as
> DecompositionPower[x,1/3] For
> > > the given function g then we can try
> > >
> > >
> DecompositionPower[x_,p_]:=Nest[ginv,Power[Nest[g,x,9],Power[2,p]],9
> > > ]
> > >
> > > and we might for instance plot
> > >
> > > Plot[DecompositionPower[0.,p],{p,0,1}]
> > >
> > > The function DecompositionPower[x_,p_] should be quite easy to
> > > generalize to a large set of polynomial functions g.
> Maybe it then
> > > should have the name FractionalNest?
> > >
> > > Best regards
> > >
> > > Ingolf Dahl
> > >
> > > -----Original Message-----
> > > From: Kelly Jones [mailto:kelly.terry.jones at gmail.com]
> > > Sent: 28 November 2006 12:04
> > > To: mathgroup at smc.vnet.net
> > > Subject: [mg71839] [mg71764] Functional decomposition (solving
> > > f[f[x]] = g[x] for given g)
> > >
> > > Given a function g[x], is there a standard methodology to find a
> > > function f[x] such that:
> > >
> > > f[f[x]] == g[x]
> > >
> > > A "simple" example (that doesn't appear to have a closed form
> > > solution):
> > >
> > > f[f[x]] == x^2+1
> > >
> > > There are probably several (approximable?
> > > power-series-expressible?) answers, but the most natural answer
> > > would be everywhere positive and monotone increasing for x>0.
> > >
> > > For x^2, there are at least two continuous solutions:
> > > x^Sqrt[2] and x^-Sqrt[2], the former being more "natural".
> > > It's somewhat amazing that x^2+1 is so much harder.
> > >
> > > Of course, the next step would be to find f[x] such that:
> > >
> > > f[f[f[x]]] == g[x]
> > >
> > > for given g[x], and so on.
> > >
> > > Thought: is there any operator that converts functional
> composition
> > > to multiplication or something similar? That would reduce this
> > > problem to find nth roots and applying the inverse operator?
> > >
> > > Other thought: For some reason, taking the geometric
> average of the
> > > identity function and g, and then iterating, seems like a
> good idea.
> > > EG, Sqrt[x*g[x]], Sqrt[x*Sqrt[x*g[x]]], and so on.
> > >
> > > --
> > > We're just a Bunch Of Regular Guys, a collective group
> that's trying
> > > to understand and assimilate technology. We feel that
> resistance to
> > > new ideas and technology is unwise and ultimately futile.
> > >
> > >
> > >
> >
> >
> >
>
- Follow-Ups:
- Re: RE: Re: RE: RE: Functional decomposition (solving f[f[x]] = g[x] for given g)
- From: Andrzej Kozlowski <akoz@mimuw.edu.pl>
- Re: RE: Re: RE: RE: Functional decomposition (solving f[f[x]] = g[x] for given g)