FuncFit and FindRoots

Is ist possible to use FuncFit on a function which requires a Findroots call to calculate a parameter which depends also on the fit coefficients?

The function looks like this:

Function EOS_Fit(w,T,P) : FitFunc
    Wave w
    Variable T
    Variable P

    //CurveFitDialog/ These comments were created by the Curve Fitting dialog. Altering them will
    //CurveFitDialog/ make the function less convenient to work with in the Curve Fitting dialog.
    //CurveFitDialog/ Independent Variables 2
    //CurveFitDialog/ T
    //CurveFitDialog/ P
    //CurveFitDialog/ Coefficients 7
    //CurveFitDialog/ w[0] = K
    //CurveFitDialog/ w[1] = Kprime
    //CurveFitDialog/ w[2] = y0
    //CurveFitDialog/ w[3] = yprime
    //CurveFitDialog/ w[4] = Cv
    //CurveFitDialog/ w[5] = S0
    //CurveFitDialog/ w[6] = F0
    wave ww             //wave ww resides in root and is labeled, which makes equation writing a lot easier for me ... also used to pass fit coefficients to FindRoots
    ww[%K] = w[0]
    ww[%Kprime] = w[1]
    ww[%y0] =w[2]
    ww[%yprime] = w[3]
    ww[%Cv] = w[4]
    ww[%S0] = w[5]
    ww[%F0] = w[6]  
    Variable V_TP
    ww[%T] = T
    ww[%P] = P
    FindRoots/Q /L=0.5 /H=3.5 V_Liq, ww
    V_TP = V_Root                                     //V_TP needs numerical solving
    variable f, a1, y, Ps, Pth, Fc, Fth, F_VT, G_TP, V_V0
    f = 1/2* ( (ww[%V0]/V_TP)^(2/3)-1 )    
    a1 = 3/2* (ww[%Kprime]-4)              
    V_V0 = V_TP/ww[%V0]                
    y = ww[%y0] + ww[%yprime]* (V_V0-1)  
    Fc = 9 * ww[%K] * 1e4 * ww[%V0] *(1/2 * f^2 + 1/3*a1*f^3)
    Fth = - ww[%S0] * (ww[%T] -ww[%T0]) - ww[%Cv] * (ww[%T] * ln(ww[%T]/ww[%T0]) - (ww[%T] -ww[%T0])) - ww[%Cv]*(ww[%T] -ww[%T0]) * (( ww[%y0]-ww[%yprime]) * ln(V_V0) + y - ww[%y0] )
    F_VT = ww[%F0] + Fc + Fth
    G_TP = F_VT + P*1e4 * V_TP 
    return G_TP / ww[%FU]

function V_Liq(ww, V_TP)  
    wave ww
    variable V_TP
    variable f, a1, y, Ps, Pth, V_V0   
    f = 1/2* ( (ww[%V0]/V_TP)^(2/3)-1 )
    a1 = 3/2* (ww[%Kprime]-4)
    V_V0 = V_TP/ww[%V0]
    y = ww[%y0] + ww[%yprime]* (V_V0-1)
    Ps = 3*ww[%K] * f * (2*f + 1)^(5/2) * (1+f*a1) 
    Pth = y/V_TP * ww[%Cv] * (ww[%T] - ww[%T0]) * 1e-4 
    return Ps + Pth - ww[%P]

I only get a "singular matrix error" or coefficients run into the constraint limits (although a solution should be possible within those limits), but I don't know if I can expect it to work at all, or is my set-up is flawed!?
Yes, you can use FindRoots in a fitting function, but you have to be careful. The solution from FindRoots has limited precision, and may not give absolutely smooth output. Consequently, you need to use an epsilon wave. First I would run your fit function outside of a curve fit to verify that it actually computes reasonable values.

Here is a description of the epsilon wave that I wrote for another customer:

Curve fitting uses partial derivatives of your fit function with respect to the fit coefficients in order to find the gradient of the chi-square surface. It then solves a linearized estimate of the chi-square surface to find the next estimate of the solution (the minimum in the chi-square surface).

Since Igor doesn't know anything about your fit function, it (he?) must use a numerical approximation of the derivatives, which are calculated by finite differences. That is, a model value is calculated at the present estimate of the fit coefficients, then each coefficient is perturbed by a small amount and the derivative is calculated from the difference.

The epsilon wave sets the size of the perturbation used to calculate the derivative. If you don't use an epsilon wave, Igor sets epsilon to:

if your coefficient wave is single precision,
eps[i] = 1e-4*coef[i], unless coef[i] is zero, in which case eps[i] = 1e-4
if your coefficient wave is double precision,
eps[i] = 1e-10*coef[i], unless coef[i] is zero, in which case eps[i] = 1e-10

There are a couple of reasons you might need to explicitly set the epsilon. One is if your fit function is insensitive to a coefficient. That is, perturbing the coefficient makes a very small change in the model value. You can get a situation where the dependence of the model is so small that floating-point truncation results in no change in the model. You have to set epsilon to a sufficiently large value that the model actually changes when the perturbation is applied. But this most likely not your problem.

Another case where you need to set epsilon is where the model is discrete or noisy. This can happen if the model involves a table look-up or a series solution of some sort. In the case of table look-up, epsilon needs to be large enough to make sure you get two different values out of the table.

In the case of a series solution, you have to stop summing terms in the series at some point. If the truncation of the series results in less than full floating-point resolution of the series, you need to make sure epsilon is large enough that the change in the model is larger than the resolution of the series. A series might include something like a numerical solution of an ODE (using IntegrateODE). It could also involve FindRoots or Optimize, each of which gives you an approximate result, and run faster if you don't demand high precision. Now the epsilon value must be set large enough that you get model values sufficient different that they are larger than the noise introduced by FindRoots.

In the case of FindRoots, you may also improve (but maybe not entirely solve) the problem by using /T to increase the accuracy of the result. That would be at the expense of longer computation time.

John Weeks
WaveMetrics, Inc.