Here is the function I want to invert:
rho(t, p) =
+ 998.20123 * p ** 0
- 192.9769495 * p ** 1
+ 389.1238958 * p ** 2
- 1668.103923 * p ** 3
+ 13522.15441 * p ** 4
- 88292.78388 * p ** 5
+ 306287.4042 * p ** 6
- 613838.1234 * p ** 7
+ 747017.2998 * p ** 8
- 547846.1354 * p ** 9
+ 223446.0334 * p ** 10
- 39032.85426 * p ** 11
- 0.20618513 * (t - 20) ** 1
- 0.0052682542 * (t - 20) ** 2
+ 0.000036130013 * (t - 20) ** 3
- 0.00000038957702 * (t - 20) ** 4
+ 0.000000007169354 * (t - 20) ** 5
- 0.000000000099739231 * (t - 20) ** 6
+ 0.1693443461530087 * p ** 1 * (t - 20) ** 1
- 10.46914743455169 * p ** 2 * (t - 20) ** 1
+ 71.96353469546523 * p ** 3 * (t - 20) ** 1
- 704.7478054272792 * p ** 4 * (t - 20) ** 1
+ 3924.090430035045 * p ** 5 * (t - 20) ** 1
- 12101.64659068747 * p ** 6 * (t - 20) ** 1
+ 22486.46550400788 * p ** 7 * (t - 20) ** 1
- 26055.62982188164 * p ** 8 * (t - 20) ** 1
+ 18523.73922069467 * p ** 9 * (t - 20) ** 1
- 7420.201433430137 * p ** 10 * (t - 20) ** 1
+ 1285.617841998974 * p ** 11 * (t - 20) ** 1
- 0.0119301300505701 * p ** 1 * (t - 20) ** 2
+ 0.2517399633803461 * p ** 2 * (t - 20) ** 2
- 2.170575700536993 * p ** 3 * (t - 20) ** 2
+ 13.53034988843029 * p ** 4 * (t - 20) ** 2
- 50.29988758547014 * p ** 5 * (t - 20) ** 2
+ 109.635566657757 * p ** 6 * (t - 20) ** 2
- 142.2753946421155 * p ** 7 * (t - 20) ** 2
+ 108.043594285623 * p ** 8 * (t - 20) ** 2
- 44.14153236817392 * p ** 9 * (t - 20) ** 2
+ 7.442971530188783 * p ** 10 * (t - 20) ** 2
- 0.0006802995733503803 * p ** 1 * (t - 20) ** 3
+ 0.01876837790289664 * p ** 2 * (t - 20) ** 3
- 0.2002561813734156 * p ** 3 * (t - 20) ** 3
+ 1.02299296671922 * p ** 4 * (t - 20) ** 3
- 2.895696483903638 * p ** 5 * (t - 20) ** 3
+ 4.810060584300675 * p ** 6 * (t - 20) ** 3
- 4.672147440794683 * p ** 7 * (t - 20) ** 3
+ 2.458043105903461 * p ** 8 * (t - 20) ** 3
- 0.5411227621436812 * p ** 9 * (t - 20) ** 3
+ 0.000004075376675622027 * p ** 1 * (t - 20) ** 4
- 0.00000876305857347111 * p ** 2 * (t - 20) ** 4
+ 0.000006515031360099368 * p ** 3 * (t - 20) ** 4
- 0.00000151578483698721 * p ** 4 * (t - 20) ** 4
- 0.00000002788074354782409 * p ** 1 * (t - 20) ** 5
+ 0.00000001345612883493354 * p ** 2 * (t - 20) ** 5
I usually use Python for my problems, so I tried sympy.solve, but it says that "No algorithms are implemented to solve the equation. "
Are you aware of a tool that could do it?
Or do you even know if there's a way to check if this function is invertible at all?
This function is from https://www.oiml.org/en/publications/bulletin/pdf/oiml_bulletin_july_2015.pdf, page 11 (numbered 9 in the doc)
I want to be able to invert it with respect to either variable.
The rho function is monotonic for a fixed t :
But not for a fixed p:
For most of my use cases, the problem is much simpler since t is always 20.
The polynom to inverse then becomes
rho(p) =
+ 998.20123 * p ** 0
- 192.9769495 * p ** 1
+ 389.1238958 * p ** 2
- 1668.103923 * p ** 3
+ 13522.15441 * p ** 4
- 88292.78388 * p ** 5
+ 306287.4042 * p ** 6
- 613838.1234 * p ** 7
+ 747017.2998 * p ** 8
- 547846.1354 * p ** 9
+ 223446.0334 * p ** 10
- 39032.85426 * p ** 11
with p in [0; 1]. This is strictly increasing on the intervalle, so it must have a solution for the inverse?
sympy.solve returns the following roots :
[
1.30502635519079,
-0.282825358240998 - 0.174919557643158*I,
-0.282825358240998 + 0.174919557643158*I,
-0.0229786659460121 - 0.441023703545107*I,
-0.0229786659460121 + 0.441023703545107*I,
0.434132824352787 - 0.596683531853141*I,
0.434132824352787 + 0.596683531853141*I,
0.877616891070605 - 0.549674307713161*I,
0.877616891070605 + 0.549674307713161*I,
1.20382263269017 - 0.326905440295712*I,
1.20382263269017 + 0.326905440295712*I
]
of which only the first one is real. What do i do with it?
Best Answer
In general, it is impossible to write down a formula for the inverse function of a polynomial with degree $\ge 5$. This result is known as Abel's impossibility theorem.
For example, if you want to solve $\rho(t,p) = y$ for $t$, there is no formula we can write down in terms of radicals that works. For this reason, you will need to resort to numerical methods, such as Newton's method, to find a solution.
It also can be difficult to prove that solutions exist. Because $\rho(t,p)$ is degree $6$ in $t$, and the $t^6$ term has a negative coefficient, we can tell that for each $p$ fixed, $\rho(t,p)$ has a maximum value, and therefore there exists $y_p\in\mathbb{R}$ such that $\rho(t,p) = y$ has a solution (not necessarily unique) if and only if $y\le y_p$. Finding $y_p$ is a simple optimization problem, so it suffices to look for solutions to $\frac{d}{dt} \rho(t,p) = 0$, which amounts to solving a degree $4$ polynomial in $t$. This is theoretically possible to do with radicals, but that's actually less efficient than using Newton's method to solve. For fixed $t$, $\rho(t,p)$ is a degree $11$ polynomial, so $\rho(t,p)=y$ always has at least one (and at most $11$) solutions for $p$.
When $t=20$, as you observe $\rho$ becomes a degree $11$ polynomial. Since it is decreasing on $[0,1]$, it has an inverse function, hence for each $y$ in its range, there is a unique $p\in[0,1]$ that solves $y=\rho(p)$. The roots don't help you too much to find an inverse function; the real root is the unique real solution to $\rho(p)=0$. Abel's theorem applies here too, so there is no closed for expression for the inverse function of $\rho$ in this case either, so you will need to resort to numeric methods here too.