1

$\def\P{\operatorname P}$

The goal of the Fourier Legendre series is to find an exact explicit series solution to $x^x=x+1$, not just $x=a+b+c+\dots$ with unknown series coefficients. Lagrange reversion was attempted, but it failed to converge each time or the series coefficients were not able to be evaluated.

To find the Fourier Legendre $\P_n(x)$ series solution of $\frac{x^x}{x+1}=1$, we transform $\frac{x^x}{x+1}$ so it contains the points $(0,0)$ and $(1,1)$:

$$f(x)=\frac65\frac{(x+1)^{x+1}}{x+2}-\frac35$$ The solution of $\frac{x^x}{x+1}=1$ is the solution of $f(x-1)=\frac35$, so $x=1+f^{-1}\left(\frac35\right)\approx1.776775$. Now we find the Fourier Legendre expansion of $\begin{cases}f^{-1}(x)&0\le x\le 1\\0&-1\le x<0\end{cases}$ to avoid any singularities, shown below:

The Fourier Legendre series usually has integration bounds of $\int_{-1}^1$, but we have integration bounds of $\int_0^1$ as we graphed $y=0$ for $-1\le x<0$:

$$f^{-1}(x)=\sum_{n=0}^\infty \left(n+\frac12\right)\P_n(x)\int_0^1 \P_n(t) f^{-1}(t)dt= \sum_{n=0}^\infty \left(n+\frac12\right)\P_n(x) a_n$$

Substituting $t\to f^{-1}(t)$ gives a Stieltjes integral:

$$a_n=\int_0^1 t\P_n\left(\frac65\frac{(t+1)^{t+1}}{t+2}-\frac35\right)\ d\left(\frac65\frac{(t+1)^{t+1}}{t+2}-\frac35\right)$$

Now expand $\P_n\left(\frac65\frac{(x+1)^{x+1}}{x+2}-\frac35\right)=\P_n\left(z-\frac35\right)$ as a Taylor series at $z_0=-\frac35$ to have a polynomial in $\frac65\frac{(t+1)^{t+1}}{t+2}$ via Legendre $\P_n^k(x)$ expansion:

$$\P_n(z)=\sum_{k=0}^n\frac{(-1)^k \P_n^k(z_0)}{\left(1-z_0^2\right)^\frac2kk!}(z-z_0)^k$$

substitute $t\to t+1$ to get integration bounds of $\int_1^2$, and factor out the $\frac 65$ from $\frac 65\frac{t^t}{t+1}$:

$$a_n=\sum_{k=0}^n\frac{(-1)^k\P_n^k\left(-\frac35\right)}{\left(1-\left(-\frac35\right)^2\right)^\frac k2k!}\left(\frac 65\right)^{k+1}\color{blue}{\int_1^2(t-1)\left(\frac{t^t}{t+1}\right)^kd \left(\frac{t^t}{t+1}\right)}\tag1$$

integrating by parts:

$$\color{blue}{\int_1^2(t-1)\left(\frac{t^t}{t+1}\right)^kd \left(\frac{t^t}{t+1}\right)}=\frac{\left(\frac43\right)^{k+1}}{k+1}-\frac1{k+1}\int_1^2\left(\frac{t^t}{t+1}\right)^{k+1}dt\tag2$$

$\int_1^2\left(\frac{t^t}{t+1}\right)^{k+1}dt$ requires another series and special functions, but that is not our focus for now. The problem is that multiple softwares fail past about the $26$th partial sum

Wolfram Alpha

At $n=26$, Wolfram Alpha starts plotting the Fourier Legendre series at stage $(1)$ like so:

integration by parts in the series coefficients, from $(2)$, slightly reduces this effect:

shown at $n=26$. However, plotting at $n=25$ of the same series gives a smooth curve as expected:

implying there is a calculation error starting around $n=26$. Maybe it is because $\int_0^1 \left(\frac{t^t}{t+1}\right)^{k+1}dt$, from $(2)$, grows into the ten thousands for $k\approx 26$ and the plotting discrepancy appears. Conversely, the series itself may diverge.

Mathematica (Wolfram Cloud)

This section uses the expansion from stage $(1)$. Unfortunately, past $n=10$ terms, computation time runs out:

enter image description here

and when $x=\frac35$, the Fourier Legendre series should give $x-1=0.776775\dots$, but it converges incorrectly past $n=26$ just like Wolfram Alpha:

$$\left[\begin{matrix}n&\text{partial sum}\\24&0.771697\\25&0.780971\\26&0.780639\\27&1.26721\\28&-0.648316\\30&6.7218\end{matrix}\right]$$

Mathematica notebook

Typing in Sum[(n+1/2) LegendreP[n,3/5] NIntegrate[t LegendreP[n,(6/5 (t+1)^(t+1)/(t+2)-3/5)] D[6/5 (t+1)^(t+1)/(t+2),t],{t,0,1}],{n,0,25}] has convergence errors:

  • “Nintegrate: Nintegrate failed to converge to prescribed accuracy after 9 recursive bisections in t near {t} = <0.929672}. Nintegrate obtained 0.000539183 and 1.01235×10 integral and error estimates”

  • “Nintegrate: Numerical integration converging too slowly: suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small”

The gradual convergence for $0\le n\le 26$ is expected, but the issue for $n>26$ is not. The series seems to diverge due to software errors, but does it really diverge?

Тyma Gaidash
  • 13,576
  • If this series works, we can solve $x^x=x+1$ exactly albeit in a slightly complicated form. – Тyma Gaidash Jun 17 '23 at 21:30
  • 1
    Have you tried computing this with other software? Have you tried computing with greater decimal precision? – Somos Jun 23 '23 at 00:56
  • You can use create a Wolfram Cloud account and use Mathematica online. You would be wise to use PARI/GP to have an alternative for computation to compare. Both software have the ability to plot functions and use arbitrary precision which is almost certainly what you need to use in this situation. – Somos Jun 23 '23 at 01:12
  • You need to make clearer what it is you are trying to do mathematically. I am confused how using Fourier-Legendre solves $\frac{x^x}{x+1}=1.$ Please give more details about what you are trying to do here. Also why the $P_n^m(x)$ is needed. – Somos Jun 24 '23 at 10:20
  • You may get faster convergence if you expand $f^{-1}((x+1)/2)$ instead so there's no sharp corner in the domain. – eyeballfrog Jun 24 '23 at 12:47
  • You already have $x\approx 1.77677504$ for which you can get arbitrary close approximation using Mathematica FindRoot[ x^x - x -1,{ x, 1.7}, WorkingPrecision->50] . Why Fourier-Legendre? – Somos Jun 24 '23 at 16:19
  • Yes, there are many numerical techniques, even ones that give a series. For example, Let $f(x) = (2+x)^{2+x}-(4+x)$, Then $f(-2+1.77677...) = -1.$ Expand $f^{-1}(x) = x/(3+\log(16)) + O(x^2)$ as a series and find $ f^{-1}(-1).$ – Somos Jun 24 '23 at 17:53
  • @Somos Your series expansion has no explicitly defined series coefficients, but instead uses big $O$ notation. This post may give an explicit series solution, explicit series coefficients etc unlike the form you wrote. – Тyma Gaidash Jun 24 '23 at 18:04
  • You can use the Lagrange inversion formula to get explicit coefficients. For an example, the Wikipedia article finds the solution of $x^p-x+z = 0$ as an explicit series in $z$. – Somos Jun 24 '23 at 18:29
  • @Somos I tried Lagrange reversion/inversion multiple times, but it diverged or went to $0$ each time. For example, try $-1+\sum\limits_{n=1}^\infty\frac1{n!}\left.\frac{d^{n-1}}{dx^{n-1}}x^{nx}\right|{x=-1}$. Also, series coefficients like $\lim{x\to1}\frac{d^{n-1}} {dx^{n-1}}\left(\frac{x-1}{x^x+x-2}\right)^n$ are very complicated. The next best technique to get a series expansion is the Fourier Legendre series. Thanks for the tip – Тyma Gaidash Jun 24 '23 at 18:43

1 Answers1

3

As I understand it, the problem is to solve $\, x^x = x+1 \,$ as an explicit series. The question uses a Fourier-Legendre polynomial approximation method. I suggest another method. Define function $$ f(x) := x^x - x - 1.$$ Notice that $\,f(1) = -1\,$ and $\,f(2) = 1.\,$ The root of the equation is $\,x_0 \approx 1.776775\,$ where $\,f(x_0) = 0.\,$ Use Lagrange inversion to get $\,x_0 = f^{-1}(0).\,$ Invert the power series for $\,f(x)\,$ centered at $\,x=2\,$ to get

$$ f^{-1}(1\!-\!x)\!=\!2\!-\!\frac{1}{3\!+\!4t}x -\frac{3+4t+2t^2}{(3+4t)^3}x^2 -\frac{81+108t+228t^2+132t^3+32t^4}{6(3+4t)^5}x^3 - \dots $$

where $\,t:=\log(2).\,$ Apply this result for $\,x=1\,$ to get

$$ x_0 = f^{-1}(0) = 2 - \frac{1}{3+4t} -\frac{3+4t+2t^2}{(3+4t)^3} -\frac{81+108t+228t^2+132t^3+32t^4}{6(3+4t)^5} - \dots $$

The rate of linear convergence is $\,1/2$. The corresponding Mathematica code is

ClearAll[DP, M, f, x, x0, xn]; DP = 20; (* decimal digit precision *)
f[x_] := x^x - x - 1;
x0 = x/.FindRoot[f[x], {x, 2}, WorkingPrecision->DP][[1]];
dn[n_] := xn[n] - x0;
(* xn[M] = f^(-1)(0) with f^(-1)(x) series truncated *)
xn[M_Integer] := (xn[M] = Module[{pow, fx, t, iser, equ, un, sols, xt, gcd},
pow = (3 + 4*t)^(2*M-1);
fx = Series[ f[x], {x, 2, M}] /. Log[2] -> t;
iser = 2 - Sum[un[n]*(1-x)^n/(3 + 4*t)^(2*n-1), {n, M}];
equ = CoefficientList[Normal@Series[Normal@Simplify[iser/. x->fx], {x, 0, M}], x];
sols = Solve[MapThread[Equal, {Table[If[n==1,1,0], {n,0,M}], equ}], Array[un, M]];
xt = (List@@Normal@Simplify[Series[iser/.sols[[1]], {x, 1, M}]] /. x->0)*pow;
gcd = PolynomialGCD@@xt;
Expand@Total[xt/gcd]*gcd/pow /. t -> N[Log[2], DP]]) /;M>2;
x0
Table[{n, Round[dn[n]*10^9], N[dn[n]/dn[n-1], 5]}, {n, 4, 20}] //Column
(* 1.7767750400970546975 *)
(*
{
 {4, 1896983, 0.36981}
 {5, 738272, 0.38918}
 {6, 297928, 0.40355}
 {7, 123527, 0.41462}
 {8, 52303, 0.42341}
 {9, 22520, 0.43056}
 {10, 9829, 0.43649}
 {11, 4340, 0.44148}
 {12, 1934, 0.44574}
 {13, 869, 0.44942}
 {14, 393, 0.45264}
 {15, 179, 0.45547}
 {16, 82, 0.45797}
 {17, 38, 0.46021}
 {18, 17, 0.46223}
 {19, 8, 0.46405}
 {20, 4, 0.46570}
}
*)

Note that Mathematica had great difficulty with InverseSeries[] for $n=17$ or above so I had to write my own code to do that for this case.

Somos
  • 37,457
  • 3
  • 35
  • 85
  • Thanks for the help. The “$\dots$” after the “$x_0=f^{-1}(0)=$“ indicate there is no expression for the series coefficients, besides something like $\lim_{x\to a}\frac{d^{n-1}}{dx^{n-1}}\left(\frac{x-a}{f(x)-f(a)}\right)^n$ and the series reversion formula is too complicated to use. The Fourier Legendre series, if it does converge correctly, solves this issue by having coefficients that can be evaluated. This question is a bit obscure and the your method is computationally better than in the question, so you will likely get the bounty. – Тyma Gaidash Jun 25 '23 at 02:43