# FuncFit and FindRoots

ChrLie

Tue, 10/18/2011 - 10:40 am

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]

End

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]

end

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]

End

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]

end

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!?

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.

support@wavemetrics.com

October 18, 2011 at 11:38 am - Permalink

Thanks for the advice! Setting epsilon values did the job!

Cheers

C

October 19, 2011 at 02:59 am - Permalink