Re: Re: RE: Re: RE: RE: Functional decomposition (solving f[f[x]] = g[x] for given g)

*To*: mathgroup at smc.vnet.net*Subject*: [mg72009] Re: [mg71996] Re: [mg71972] RE: [mg71885] Re: [mg71872] RE: [mg71839] RE: [mg71764] Functional decomposition (solving f[f[x]] = g[x] for given g)*From*: Andrzej Kozlowski <akoz at mimuw.edu.pl>*Date*: Fri, 8 Dec 2006 06:17:57 -0500 (EST)*References*: <200612061104.GAA02807@smc.vnet.net> <B2493AAB-0463-49C0-9E90-118BB276E47F@mimuw.edu.pl> <200612071126.GAA28291@smc.vnet.net>

A little after posting that I noticed that almost the same argument provides another counter-example: any smooth real function g with g(0) = 0 and g'(0)<0 cannot be represented as f(f(x)), with f a real function of class C1. Andrzej Kozlowski On 7 Dec 2006, at 20:26, Andrzej Kozlowski wrote: > I meant "class C1", of course: differentiable, with continuous > derivative. > > Andrzej Kozlowski > > > On 7 Dec 2006, at 17:58, Andrzej Kozlowski wrote: > >> To my mind an equally interesting problem (though admittedly not if >> one wishes to apply Mathematica) is to find the widest class of >> functions g for which the equation admits no solutions. There is >> one case which is immediately obvious: the functional equation f(f >> (x))==g(x) cannot have any solutions f (of class C2) for smooth >> functions g whose derivative is g'(x) <0 everywhere; g(x) = exp(-x) >> is such a function. The proof is completely trivial, but as far as >> I can tell no-one seems to have mentioned it explicitly so here it >> is: >> if f(f(x)) = g(x) then f'(f(x)) f'(x)) = g'(x) is negative. But >> that means that f' has opposite sign at x and f(x). That in turn >> implies (if f is of class C2) that f'(t) = 0 for some t between x >> and f(x). But that would mean that g'(t) = 0, contradicting the >> fact that g'(x)<0 everywhere. >> I have given very little thought to this but it seems to me that >> the assumption g'(x)>0 for all x should suffice for the existence >> of f (though may not suffice for the existence of a nice enough f, >> e.g. analytic). I would expect also that examples can be >> constructed where g'(x)>0 for some x and g'(y)< 0 for some why and >> where no f of class C2 can exists, but I can't find any obvious ones. >> >> Andrzej Kozlowski >> >> On 6 Dec 2006, at 20:04, Ingolf Dahl wrote: >> >>> 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. >>>>>> >>>>>> >>>>>> >>>>> >>>>> >>>>> >>>> >>> >>> >> >

**References**:**RE: Re: RE: RE: Functional decomposition (solving f[f[x]] = g[x] for given g)***From:*"Ingolf Dahl" <ingolf.dahl@telia.com>

**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)**

**FullSimplify and HypergeometricPFQ**

**Re: RE: Re: RE: RE: Functional decomposition (solving f[f[x]] = g[x] for given g)**

**Re: RE: Re: RE: RE: Functional decomposition (solving f[f[x]] = g[x] for given g)**