# Solve functional equation $h(y)+h^{-1}(y)=2y+y^2$

I was trying to solve a certain physics problem, and encountered the functional equation that contains a function $$h$$ and its inverse $$h^{-1}$$: $$$$h(y)+h^{-1}(y)=2y+y^2.\tag{1}$$$$

Q: Does $$(1)$$ have a unique solution and is it possible to find it in closed form?

Equation $$(1)$$ looks quite simple and probably it is simple to analyze too, but I couldn't figured out how.

Background information: The physics problem I was trying to solve is finding the dependence of the current $$J_0$$ on the voltage $$U_0$$ in this infinite chain that contains nonlinear elements with quadratic volt-ampere characteristics $$I(V)=\alpha V^2$$ and ohmic resistors $$R$$:

According to dimensional analysis one can write $$J_0(U_0)=\frac{1}{\alpha R^2}f(\alpha RU_0).$$ Solving simple system of equations I obtained the following functional equation for the unknown function $$f$$ $$f(x)=(x-f(x))^2+f(x-f(x)).\tag{2}$$ Now introduce another function $$h$$ according to $$x-f(x)=h(x).$$ Then equation $$(2)$$ becomes $$x-h(x)=h^2(x)+h(x)-h(h(x)).$$ Let $$h(x)=y$$, then $$x=h^{-1}(y)$$ and we get $$(1)$$.

This is not a textbook problem and I don't even know whether it has a solution. I was studying it out of curiosity.

• See this answer math.stackexchange.com/q/3317234 Dec 8 '19 at 11:15
• @Nemo just read your question. can you please make it good? add domain and range of $h$. make it a proper math question. are you assuming $h$ is a bijection (so that $h^{-1}$ makes sense)? are you assuming continuity? Dec 17 '19 at 3:15
• @Nemo dude. what is the domain. where is this thing defined? im talking about $h$ only Dec 17 '19 at 9:25
• If $x=0$, then there's no external potential and one has zero current, i.e., $f(0)=0$. Negative potentials ($x<0$), on the other hand, are unphysical because the volt-ampere characteristic for the nonlinear element always gives a positive current. (By contrast, the volt-ampere characteristic $I(V)=\beta V^3$ should give a chain which is physical for both positive and negative external potentials.) So the domain is $x\geq 0$ but the case of $x=0$ is trivial. Dec 17 '19 at 14:40
• One additional reason to be interested in the case $I(V)=\beta V^3$ is that the resulting difference equation, $$u_{k+1}-2u_k+u_{k-1} = u_k^3,$$ is an instance of the (stationary) discrete nonlinear Schrodinger equation. See page 4 of macs.hw.ac.uk/~chris/em02.pdf. Dec 17 '19 at 16:52

Unfortunately, your proposed functional equation: $$h(y)+h^{-1}(y)=2y+y^2 \tag{1}$$ has no differentiable solutions near $$\,y=0.\,$$ Suppose that $$h(y) = 0 + a_1 y + O(y^2). \tag{2}$$ Then its inverse function is $$h^{-1}(y) = 0 + \frac1{a_1} y + O(y^2). \tag{3}$$ Substituting equations $$(2)$$ and $$(3)$$ into $$(1)$$ gives $$h(y)\!+\!h^{-1}(y) = \frac{a_1^2+1}{a_1}y \!+\! O(y^2) = 2y \!+\! O(y^2). \tag{4}$$ The only solutions for $$\,a_1\,$$ are $$\,a_1=1\,$$ and $$\,a_1=-1.\,$$ However, we are given that $$0 and thus the only possible value is $$\,a_1=1.\,$$ Suppose we want more terms in the power series $$h(y) = 0 + y + a_2y^2 + O(y^3). \tag{6}$$ The inverse function is now $$h^{-1}(y) = 0 + y - a_2y^2 + O(y^3). \tag{7}$$ Adding these two equations together is the equation $$h(y) + h^{-1}(y) = 0 + 2y + 0y^2 + O(y^3). \tag{8}$$ This implies that equation $$(1)$$ can not be satisfied.

There is a possibility that there exists an exponent $$\,e\,$$ not an integer such that $$\, h(y) = 0 + y + c y^e + \cdots \,$$ which perhaps should be studied according to comments.

In fact, define $$g(x) : = \sqrt{h(x^2)}. \tag{9}$$ Then solving for a power series expansion gives $$g(x) = x - \frac{x^2}{\sqrt{6}} + \frac{x^3}6 + O(x^4) \tag{10}$$ which implies that $$h(x^2) = (x^2+x^4/2)+\frac1{\sqrt{6}}f(x) \tag{11}$$ where $$f(x) \!=\! -2 x^3 \!-\!\frac{11x^5}{24} \!+\! \frac{117x^7}{1280} \!-\!\frac{5491x^9}{110592} \!+\! \frac{156538363x^{11}}{3715891200} \!+\! O(x)^{13}. \tag{12}$$

A Wolfram Language code to calculate $$\,g(x)\,$$ is

ClearAll[x, g, gx];
gx[3] = x - 1/Sqrt[6]*x^2 + O[x]^3;
Do[g = Normal[gx[n]] + O[x]^(2+n); gx[n+1] = Simplify[
g + (x^2 + (Normal[g]/.x -> g)^2 - g^4 - 2*g^2) * 3 /
((4+n)*x^2*Sqrt[6])],
{n, 3, 6}]


As some comments suggest, it seems that the power series for $$\,f(x)\,$$ has zero radius of convergence. This makes finding properties of it difficult. Perhaps we need to let $$\,x\,$$ approach $$\,\infty.$$ In this case we find that $$\,h(x) \approx \sqrt{x}\,$$ with infinite number of other terms. Using equation $$(1)$$ we can find the expansion $$h(x) = x^{1/2} - 1 + \frac12 x^{-1/4} + \frac1x s(x) \tag{13}$$ where $$s(x) := \sum_{k=2}^\infty 2^{-k}(x^{-\frac32 2^{-k}} - x^{-2^{-k}} ). \tag{14}$$ This gives moderately good approximations as $$\,x\,$$ gets large and down to $$1.$$

One method that leads to equations $$(13)$$ and $$(14)$$ is as follows. From equation $$(1)$$ we immediately get $$h^{-1}(y) = y^2 + 2y - h(y) \tag{15}$$ and if $$\,y=h(x),\,$$ then $$x = h(x)^2 + 2h(x) - h(h(x)). \tag{16}$$ We start with an approximation and try to find what additional term will satisfy equation $$(16)$$. So we guess $$h(x) = x^{1/2} + cx^e + \cdots \tag{17}$$ where $$\,\cdots\,$$ denotes terms with smaller exponents. Substitute equation $$(17)$$ in equation $$(16)$$ to get $$x \!=\! (x \!+\! 2cx^{1/2+e} \!+\! \cdots) \!+\! 2x^{1/2} \!+\! \cdots \!=\! x \!+\! 2x^{1/2}(cx^e \!+\! 1) + \cdots. \tag{18}$$ This implies $$\,0 = cx^e + 1.\,$$ Solving this we get our next guess as $$\,h(x) = x^{1/2} - 1 + cx^e + \cdots. \tag{19}$$ Repeating this process leads to equations $$(13)$$ and $$(14)$$.

The series in $$(14)$$ appears to converge but I have no proof, only numerical evidence. I also have no proof that $$(13)$$ satisfies the functional equation $$(1)$$.

An answer to this question by 'Semiclasical' contains sequence recursions $$u_n = u_{n-1} - j_{n-1} \quad \text{ and } \quad j_n = j_{n-1} - u_n^2. \tag{20}$$ with the property that $$\, u_n = h(u_{n-1})\,$$ and that $$\,u_n \to 0^+.\,$$ I found that for suitable starting values of $$\,u_1\,$$ and $$\,j_1\,$$ we have $$u_n = 6n^{-2} - \frac{15}2n^{-4} + \frac{663}{40} n^{-6} - \frac{43647}{800}n^{-8} + \cdots. \tag{21}$$

An answer to this question by 'Cesareo' uses recursion to construct a sequence of functions $$\,h_n(x)\,$$ which seems to converge to a global solution. It would be nice to give a proof of this.

• I think this dovetails with a point I made in comments earlier: the quadratic nonlinearity of the elements ensures that the model is only physical for nonnegative external potentials. In particular, this means that the domain and codomain of $h$ are the nonnegative reals. (It also fits with the computations in my answer which suggest $h(x)\sim x^{1/2}$.) Dec 18 '19 at 0:04
• One correction: The asymptotics I get are $h(x)=x-x^{3/2}/\sqrt{3}$. So while $h'$ exists for $x\geq 0$, it certainly doesn't for $x<0$. Dec 18 '19 at 0:46
• @Semiclassical Thanks for an interesting comment! Dec 18 '19 at 1:00
• @Semiclassical , Somos, is it possible to calculate more coefficients of this series? It seems that coefficients of integer powers might all be rational, and coefficients of half integer powers are $\sqrt{6}$ times rationals. If this is indeed true then one can form integer sequences from these coefficients and look them up in oeis... Dec 18 '19 at 7:42
• I was reviewing this question and noticed eqs. 13-14 that I missed. Seems like this series converges fast for even small $x$. How did you obtain this equation? Do you have a proof? Mar 30 '20 at 10:07

From

$$x = 2h(x)+h^2(x)-h(h(x))\Rightarrow h(x) = \pm\sqrt{h(h(x))+x+1}-1$$

now using an iterative approximation over the positive leaf such as

$$h_{k+1}(x) = \sqrt{h_k(h_k(x))+x+1}-1$$

or according to the MATHEMATICA script which follows, (inspired after a fruitful discussion with Semiclassical)

Clear[h]
h[x_, 1] := x
h[x_, n_] := h[x, n] = Sqrt[h[h[x, n - 1], n - 1] + x + 1] - 1


we obtain beginning with $$h_1(x) = x$$ the successive approximations in red, and black for $$h_5(x)$$ showing a good convergence for $$x > 0$$.

Also the accomplishment for $$h^{-1}(x)+h(x) = 2x+x^2$$ can be analyzed in the following plot

err5 = y - h[2 y + y^2 - h[y, 5], 5];
Plot[Abs[err5], {y, 0, 4}, PlotStyle -> {Thick, Black}]


• Would I be correct in assuming that the convergence for $x>2$ is not good? That would be a disappointment---the large-$x$ behavior is not uninteresting---but not at all a shock. Dec 17 '19 at 15:08
• Your Mathematica code is not doing what you think it is doing. In particular, the code to set h[x_] is wrong. The fourth degree polynomial should be $\,x-x^2+3/4x^3-3/16x^4\,$. Instead you have $\,x-3/4x^2+x^3/2-x^4/8\,$ which is wrong! Dec 19 '19 at 12:37
• I'll note that my own code uses a different construction: Clear[h]; h[x_, 1] := x; h[x_, n_] := h[x, n] = 1/2 (x - h[x, n - 1]^2 + h[h[x, n - 1], n - 1]);. What I obtain for $h_3(x)$ is in agreement with Somos as well as some other numerical checks. Dec 19 '19 at 19:26
• For specificity, I have the following sequence: $h_1(x)=x$, $h_2(x) = \frac12 (x-x^2+x)=x-x^2/2$, and then $$h_3(x) = \frac12\left[x-h_2(x)^2+h_2(h_2(x))\right]=\frac12 \left[x-(x-x^2/2)^2+((x-x^2/2)-(x-x^2/2)^2/2)\right].$$ Dec 19 '19 at 19:53
• You don't really need to plot anything to check that. Is the recursion I wrote correct? If so, then the highest-order terms are $$\frac12 \left [-(-x^2/2)^2-(-x^2/2)^2/2\right ] =-\frac{1}{8}x^4 -\frac{1}{16}x^4 = -\frac{3}{16}x^4\neq -\frac{x^4}{8}.$$ Dec 19 '19 at 19:59

The following is not intended as an answer so much as it is a development of the physics problem. (Hence why it's a wiki answer.) Let $$U_k$$ be the potential difference across the $$k$$th nonlinear element with corresponding current $$\alpha U_k^2$$, for $$k\geq 1$$. By virtue of Kirchoff's loop rule, we deduce that the voltage across the $$k$$th resistor is given by $$U_{k-1}-U_{k}$$. Since each resistor is ohmic, we conclude that the current in the $$k$$th resistor is $$J_{k-1}=R^{-1}(U_{k-1}-U_{k}).$$ In particular, we have $$J_0=R^{-1}(U_{0}-U_{1})$$. Kirchoff's junction rule then demands $$J_{k} = J_{k-1}-\alpha U_k^2.$$ Introducing $$(u_k,j_k):=(\alpha R U_k,\alpha R^2 J_k)$$ and rearranging, we obtain the dimensionless recurrence relations \begin{align} u_k = u_{k-1}-j_{k-1},\qquad j_k &= j_{k-1}-u_k^2\\&=j_{k-1}-(u_{k-1}-j_{k-1})^2. \end{align} This may be compactly expressed as $$(u_k,j_k)=g(u_{k-1},j_{k-1})=\cdots=g^k(u_0,j_0)$$ where $$g(u,j):=(u-j,j-(u-j)^2)$$. That is, the sequence of voltages and currents is obtained by iterating from an initial choice of $$(u_0,j_0)$$.

However, on physical grounds we are only concerned with solutions for which the currents and voltages are positive and monotonically decreasing to zero. This is rather delicate, as numerical experimentation demonstrates that the fixed point $$(u,j)=(0,0)$$ is badly unstable. (I don't know how to formally prove this. We can note, though, that the conditions $$u-j>0$$ and $$j>j-(u-j)^2>0$$ require $$j. Hence the range of possible $$j$$ gets smaller and smaller as $$j\to 0^+$$, which to me seems consistent with the origin being unstable.)

Finally, to obtain the advertised functional equation, we seek a solution of the form $$j_k=f(u_k)$$. In that case, applying the first equation to $$j_{k+1}=f(u_{k+1})$$ yields

$$j_{k+1} = f(u_{k+1}) = f(u_k-j_k) = f(u_k-f(u_k)).$$

On the other hand, from the second equation we have

$$j_{k+1} = j_k - u_{k+1}^2 = j_k - (u_k-j_k)^2 = f(u_k)-(u_k - f(u_k)^2.$$ Together, we obtain the desired functional equation $$f(x-f(x)) = f(x)-f(x-f(x))^2$$ under the identification $$x=u_k$$. (From this we may further conclude that $$h(u_k) = u_{k+1}$$, i.e., the voltages are obtained by iterating $$h$$ on $$u_0$$.)

I'm not confident I know how to justify the condition $$j_k=f(u_k)$$ rigorously. Physically, however, it has a simple interpretation: Suppose we draw a vertical line between the first nonlinear element and the second resistor. Then on the right we have a copy of the infinite chain, but subject to voltage $$V_1$$ and $$J_1$$. Since it's the same chain, it must have the same voltage-current relationship as the original, i.e., $$j_1=f(u_1)$$. We may then repeat this logic with the next element-resistor pair and so on, yielding $$j_j=f(u_k)$$ for all $$k$$ as desired.

At this point I altogether run out of firm conclusions. But I do have some more observations:

• In a problem with multiple dimensionful parameters, it's often wise to study limiting cases for which the problem simplifies. For instance, one has the trivial limit $$U_k,J_k\to 0$$ as $$U_0\to 0$$. More interesting is the case $$\alpha\to \infty$$, where each nonlinear element is a short circuit and so the entire current $$J_0$$ will flow through the first nonlinear element (path of least resistance). Therefore $$J_0=R^{-1} U_0$$ as $$\alpha R U_0\to\infty$$. Other limits are not helpful: If $$\alpha\to 0$$, then all the nonlinear elements become open circuits and therefore one has an infinite chain of identical resistors, i.e., infinite resistance. Similarly, if $$R\to 0$$ then the voltage across each nonlinear element is $$U_0$$ and therefore the current required would be infinite. Neither situation is physical and so there's no evident conclusion to draw. (To render them physical, one could either introduce ohmic resistances on the branches or consider a finite chain.)

• It's worth noting that, while the dimensionless variables chosen above seem obvious enough, they're not the only ones possible. For instance, we could just have well have taken $$u'_k:=U_k/U_0$$, $$j'_k:=(R/U_0)J_k$$ to obtain the dimensionless equations $$j'_{k-1} = u'_{k-1} - u'_k,\qquad j'_k = j'_{k-1} - \gamma {u'_k}^2$$ where $$\gamma = \alpha R U_0$$. The principal benefit of this is that we can now explicitly consider $$\alpha \to \infty$$, in which case $$\gamma\to \infty$$ and so the second equation collapses to $$u'_k=0$$ for $$k\geq 1$$. Hence $$j'_1=u_0$$ as stated previously. This moreover raises the possibility of solving the equations perturbatively in powers of $$\gamma^{-1}$$, though I've run into trouble proceeding along this line.

• In the above, I've presented the problem in terms of a system of nonlinear first-order difference equations. This is readily converted into a single nonlinear difference equation of second-order: $$u_k^2 = j_{k-1}-j_k = (u_{k-1}-u_k)-(u_{k}-u_{k+1})=u_{k+1}-2u_k+u_{k-1}.$$ This second-order difference equation is analogous to the differential equation $$u(x)^2 = 2u''(x)$$ with $$j(x):=-u'(x)$$. By taking a first integral, we obtain $$\frac{1}{3}u(x)^3 = u'(x)^2+C.$$ The requirement that $$u(x),j(x)\to 0$$ as $$x\to \infty$$ then imposes $$C=0$$. The resulting first-order ODE is separable with solution $$u(x)=u(0)(1+x\sqrt{u(0)/12})^{-2}$$. Then $$j(0)=u(0)^{3/2}/\sqrt{3}$$, suggesting $$j_0\propto \frac{u_0^{3/2}}{\sqrt{3}}\implies J_0\propto \frac{(\alpha R U_0)^{3/2}}{\sqrt{3}\alpha R^2}$$ (for small $$U_0$$, I think?) as the asymptotics for the original difference equation. The above is not exactly rigorous, so I don't know how much stock to put in it; however, it does seem to match up with what numerics I've done.

• One solution idea, in line with the above use of a first integral to solve the differential equation, is to search for a conservation law, i.e., a function $$H(x,y)$$ such that $$H(u_k,j_k) = (H\circ g)(u_{k-1},j_{k-1})=H(u_{k-1},j_{k-1}).$$ (Or, for $$u_k$$ alone, a function $$H'$$ such that $$H'(u_{k+1},u_k)=H'(u_k,u_{k-1})$$.) This would dictate $$H(u_0,j_0)=\cdots =H(u_k,j_k)=\cdots=\lim_{k\to \infty} H(u_k,j_k)=H(0,0).$$ As such, the set of physical $$(u_0,j_0)$$ would be prescribed as the level set of $$H$$ through $$(0,0)$$. However, I've yet to come up with such $$H(x,y)$$ and so have no definitive conclusions here.

• Nice analysis, thanks. Dec 16 '19 at 17:56
• @Nemo Thanks! I've also gone and added the derivation to the functional equation (which I was somehow not seeing until just now). Dec 16 '19 at 19:10
• And now I've added the physical meaning to the condition $j_k=f(u_k)$, which was previously eluding me. Dec 16 '19 at 21:26
• An additional thought that I don't know what to do with yet: The (dimensionless) power dissipated by the $k$th resistor is $j_{k-1}(u_{k-1}-u_k)$, while the power dissipated by the $k$th nonlinear element is $u_k^2(u_k)=(j_{k-1}-j_k)u_k$. Together, their combined power dissipated takes the form $j_{k-1} u_{k-1}-j_{k}u_k.$ Summing over $k$ and requiring that the total power dissipated be finite, we obtain $j_0 u_0$ as the total power. That's as expected, but my hope was that this would give insight into a conservation law... Dec 17 '19 at 5:03

I tried to find a more general chain that has an exact solution and then obtain the chain in the question as a limiting case of this more general problem. However, it turned out that for exact solution of this more general problem there needs to be a constraint for parameters. Because of this constraint this approach didn't work out.

Consider the infinite network in the picture:

Here $$I_X(V)=\alpha\sqrt{V}$$ and $$I_Y(V)=\beta V^2$$ are two nonlinear elements. One can combine the constants $$R$$, $$\alpha$$, $$\beta$$ to obtain a quantity that has dimension of voltage $$$$V_0=(\alpha/\beta)^{2/3},$$$$ and a dimensionless parameter $$$$\varepsilon=R(\alpha^2\beta)^{1/3}.$$$$ Then $$$$\alpha=\varepsilon\sqrt{V_0}/R,\quad \beta=\varepsilon/(V_0R).$$$$

To simplify calculations let's introduce shorthand $$Z$$ for the part of the network as shown in picture a) below.

One can easily calculate its equivalent volt-ampere characteristic $$$$\label{Z} I_Z(V)=\frac{V_0}{R}F_Z\left(\frac{V}{V_0}\right),\quad \text{где} \quad F_Z(t)=\varepsilon t^2+\varepsilon\sqrt{\varepsilon^2+t}-\varepsilon^2.$$$$

According to dimensional analysis we can write for the VAC of the infinite chain $$$$\label{AB} I_{AB}(U_0)=\frac{V_0}{R}F\left(\frac{U_0}{V_0}\right),$$$$ where $$F$$ - is an unknown dimensioneless function of a dimensionless argument. We want to find this function.

We find the following equation for $$F$$: $$$$\label{f} F(t)=F(t-F(t))+F_Z(t-F(t)).$$$$We want to find a solution $$0.

Introducing $$$$\label{h} H(t)=t-F(t),$$$$ we get $$$$\label{hfinal} H(t)+H^{-1}(t)=2t+F_Z(t).\tag{*}$$$$

Assume the solution is of the form $$F(t)=t+\frac{R}{V_0}(a-\sqrt{a^2+bV_0t})$$ with some constants $$a$$ and $$b$$.

Then $$H(t)=\frac{R}{V_0}(\sqrt{a^2+bV_0t}-a),$$ and $$H^{-1}(t)=\frac{V_0t^2+2aRt}{bR^2}.$$

Substituting these into $$(*)$$ we get $$\frac{R}{V_0}(\sqrt{a^2+bV_0t}-a)+\frac{V_0t^2+2aRt}{bR^2}=2t+\varepsilon t^2+\varepsilon\sqrt{\varepsilon^2+t}-\varepsilon^2.\tag{**}$$

One can easily see that in order for this equation to satisfy identically for all $$t>0$$ one needs to satisfy the conditions $$\varepsilon=1,\quad a=V_0/R, \quad b=V_0/R^2.$$

Unfortunately, replacing the resistance $$2R$$ with more general $$r$$ doesn't lead to anything new, because the requirement analogous to $$(**)$$ can only be satisfied when $$r=2R$$.

Thus we get that we can find the VAC of this infinite chain under the constraint $$\alpha^2\beta R^3=1$$: $$I_{AB}(U_0)=\frac{U_0+V_0-\sqrt{V_0^2+V_0U_0}}{R}, \quad \text{where} \quad V_0=(\alpha/\beta)^{2/3}.$$

Unfortunately, the requirement $$\alpha^2\beta R^3=1$$ spoils everything and it is not possible to use this solution to find the solution of the initial question as a limiting case.