I have two sets of experimental data X and Y. The data in both sets are related by $y = f(x),\hspace{0.2cm} x \in X,\hspace{0.2cm}y \in Y$ the function $f(x)$ is unknown. The data is in the interval $x,y\in[0,1]$ or better put $X\in[0,1]$ and $Y\in[0,1]$. The relation between $X$ and $Y$ is highly non linear, and by highly nonlinear (I know that is not a mathematical term) I mean that the function can't be properly fitted by polynomial series like Taylor series, Legendre polynomials, Laguerre polynomial, Hermite polynomials, I even try Bernoulli polynomials. The function apparently can't be fitted by Fourier series either. I'm currently trying wavelets techniques but as far as I get the results apparently say that it can't be fitted by wavelets either. The function $f(x)$ although closed in the interval $[0,1]$ shows a extremely non linear behavior given that the derivatives approximations are high valued and change in sign in some points. The function don't oscillate much, in some sub intervals looks like a line. I have to say that the quality of data is very good. There is no doubt on the values of the data.
The values of the underlying function can't just be approximated because the values have a precisely physical meaning and a significance, that means that small error in the approximations incurred in a big error in the physical model. That's the reason why I need not a good, but a very good approximation to the function values. Also I need that the resulting fitted model can be reasonable computed in a normal desktop machine. After doing the fitting the resulting function is going to be massively evaluated in another algorithm so I need the most efficient evaluation of $f(x)$.
The question finally is: there is another advanced or new technique to perform a curve fitting to data with these characteristics?
I do not think that the fitting model would be simple. We have $$y=\frac{V_s-V_t}{V_c-V_t} \qquad \text{and} \qquad x=\frac{P_s-P_t}{P_c-P_t}$$
Assume that we use the simplest possible models (using $T_r=\frac T{T_c}$).
and we know the values of $(V_c,V_t,P_c,P_t)$. Introducing reduced volume and pressure at saturation $$V_r=\frac {V_s}{V_c} \qquad \text{and} \qquad P_r=\frac {P_s}{P_c}$$
$$y=\frac{V_r-a}{1-a} \qquad \text{and} \qquad x=\frac{P_r-b}{1-b}$$ where $a=\frac{V_t}{V_c}$ and $b=\frac{P_t}{P_c}$ are known constants.
So, from $V_r$, we have $$T_r=1-\frac{\log ^3(V_r)}{\log ^3(Z_c)}$$ and, from $P_r$ $$T_r=\frac{k}{k-\log (P_r)}$$
Making them equal would lead to $$\log(V_r)=-\sqrt[3]{\frac{\log (P)}{k-\log (P)}} \log(Z_c)$$
In practice, I think that is would be much safer to fit independently $V_s$ and $P_s$ as function of $T$.
Edit
Reworking what we did at the time of Nist/Asme project about water properties, I reworked the problem as stated in the question. What I did is the following $$y=\frac{V_s-V_t}{V_c-V_t}=\frac{1}{1+\sum _{n=1}^p a_n\, \tau ^{n/3}}\qquad \text{where} \qquad\tau=1-T_r$$ Using $p=6$ gives an $R^2=0.999879$ and a maximum absolute error equal to $0.004$. Notice that this could be improved a lot (using a stepwise regression for finding the best exponents).
As a function of $\tau$ the function $y$ is quite nice. So, what I suggest is to use as a model $$y=\frac{1}{1+\sum _{n=1}^p a_n\, \Big[\tau(x)\Big] ^{n/3}}$$ and it is quit simple to make a model (inversion of the vapor pressure) $$\tau(x)= F[x]$$