How to approximate an implicit function using simple functions (e.g. polynomials)?

Question

I am trying to approximate the following function using taylor series (to express $y$ as an explicit function of $x$ for $0.3<y<1$):

$$ x = a (1-y)^b + c \sinh(d (1-y)) $$

where $a$, $b$, $c$, and $d$ are constants. The following is my code:

ClearAll["Global`*"]

(*Variables*)
a = 268877.11259000533`;
b = 33;
c = 0.8333333333333333`;
d = 1.5` ;
x[y_] = a*(1 - y)^b + c*Sinh[d*(1 - y)];

(*Visualization*)
Plot[x[y], {y, 0.3, 1}, PlotRange -> {0, All}, 
 PlotLabel -> "x as a function of y"]

(*Approximating using taylor series*)
taylorY[x_] = 
  Normal[InverseSeries[Series[x[y], {y, 1, 50}]]] /. {y -> x};
Plot[taylorY[x], {x, 0, 3}, PlotRange -> {0.3, 1}, 
 PlotLabel -> "y as a function of x (Taylor)"]

As shown above, the results do not converge beyond $x > 0.7$. I have tried centering the taylor series about different $y$ values, and it still couldn't fully converge within my range of interest ($0.3<y<1$).

Is there idea to fix this issue?

Answer 1

I made a WFR function to do this called PolynomialApproximation:

ClearAll["Global`*"]
(*Variables*)
a = 268877.11259000533`;
b = 33;
c = 0.8333333333333333`;
d = 1.5`;
x[y_] := a*(1 - y)^b + c*Sinh[d*(1 - y)];

approx = 
 ResourceFunction["PolynomialApproximation"][
  a*(1 - y)^b + c*Sinh[d*(1 - y)], {y, 0.3, 1}, 20, "Chebyshev"]

(*Visualization*)
Plot[{x[y], approx}, {y, 0.3, 1}, PlotRange -> {0, All}, 
 PlotLabel -> "x as a function of y and approximation"]
Plot[x[y] - approx, {y, 0.3, 1}, PlotRange -> All, 
 PlotLabel -> "difference"]

giving:

Answer 2

$Version

(* "14.3.0 for Mac OS X ARM (64-bit) (July 8, 2025)" *)

ClearAll["Global`*"]

(*Variables*)
a = 268877.11259000533`;
b = 33;
c = 0.8333333333333333`;
d = 1.5`;
x[y_] = a*(1 - y)^b + c*Sinh[d*(1 - y)];

plt = ParametricPlot[{x[y], y}, {y, 3/10, 1},
   PlotRange -> All];

data = Cases[plt, Line[pts_] :> pts, All][[1]];

y[x_] = FindFormula[data, x]

(* 0.996004 - 0.559313 x - 1.18396 x^2 + 2.01856 x^3 - 1.17682 x^4 + 
 0.304992 x^5 - 0.0296955 x^6 *)

Legended[
 Show[plt,
  Plot[y[x], {x, 0, 3.1}, PlotStyle -> {ColorData[97][2], Dashed}],
  AxesLabel -> (Style[#, 14] & /@ {x, y}),
  AspectRatio -> 1/2],
 Placed[
  LineLegend[ColorData[97] /@ {1, 2}, {"Implicit", "Explicit"}],
  {.5, .5}]]

Answer 3

Least squares is probably a decent choice:

data = Table[{y, x[y]}, {y, 0.3, 1, .01}]

polynomial (11 terms is probably more than enough)

model = Sum[coeff[i] y^i, {i, 0, 10}]

the coefficients to be obtained

coeffs = coeff /@ Range[0, 10]

least squares fit:

fit = FindFit[data, model, coeffs, y]

the approximating polynomial:

fittedModel = model /. fit

visual check:

Show[
 ListPlot[data, PlotStyle -> Red],
 Plot[fittedModel, {y, .3, 1}]
 ]

Answer 4

A continued fraction expansion also seems like it can do a decent job here.

a = 268877.11259000533`;
b = 33;
c = 0.8333333333333333`;
d = 1.5`;
F[y_, a_, b_, c_, d_] := a*(1 - y)^b + c*Sinh[d*(1 - y)];

First generate some data which we will use to fit the model.

y = Table[i, {i, 0.3, 1, .01}]
x = F[y, a, b, c, d]

Now we can fit the model.

fittedModel = 
 NonlinearModelFit[Transpose[{x, y}], 
  a1 - a2* x1 + a3 *x1^2 - a4 x1^4 + a5/(x1 + 1/(a6 + x1)) - a7/(
   x1^4 + 1/(a8 + x1^4)) + a8/(
   x1^2 + 1/(a9* x1^2 + 1/(a10 + x1^4))), {a1, a2, a3, a4, a5, a6, a7,
    a8, a9, a10}, x1, MaxIterations -> 10000, Method -> "NMinimize"]

This results in the following model.

Normal[fittedModel]

$1.32564 - 0.697703 x1 + 0.164427 x1^2 - 0.00328264 x1^4 - 0.402326/( x1 + 1/(0.806882 + x1)) - 0.11395/( x1^4 + 1/(-0.173279 + x1^4)) - 0.173279/( x1^2 + 1/(2.01633 x1^2 + 1/(8.78212 + x1^4)))$

Plotting the original data vs the fitted model we see

Show[Plot[fittedModel[x1], {x1, 0, 3}, PlotStyle -> Red], 
 ListPlot[Transpose[{x, y}], PlotStyle -> {PointSize[0.012]}], 
 Frame -> True, FrameLabel -> {"X", "Y"}]