MATLAB Answers

Moj
0

nonlinear curve fitting a function on data

Asked by Moj
on 9 Jun 2019
Latest activity Edited by Matt J
on 10 Jun 2019
Hi,
I have some experiment data. Hereby, I need to fit the following function to determine one of the variable. A Levenberg–Marquardt least-squares algorithm was used in this procedure.
My experiment data:
beta = 1.135e-4;
sin(theta) = [-0.81704 -0.67649 -0.83137 -0.73468 -0.66744 -0.43602 0.45368 0.75802 0.96705 0.99717 ];
x = [72.01 59.99 51.13 45.53 36.15 31.66 30.16 29.01 25.62 23.47 ];
And function which is needing to be fitted:
sin(theta) = -1+2*sqrt(alpha/x)*exp(-beta*(x-alpha)^2)
Is there any suggestion to determine alpha variable?
Thanks in advance

  0 Comments

Sign in to comment.

2 Answers

Answer by John D'Errico
on 9 Jun 2019
Edited by John D'Errico
on 9 Jun 2019
 Accepted Answer

sin(theta) is not a variable name.
beta = 1.135e-4;
sintheta = [-0.81704 -0.67649 -0.83137 -0.73468 -0.66744 -0.43602 0.45368 0.75802 0.96705 0.99717 ];
x = [72.01 59.99 51.13 45.53 36.15 31.66 30.16 29.01 25.62 23.47 ];
Before you do ANYTHING, PLOT your data. Ask, if it is consistent with your model. The answer here, is not really.
plot(x,sintheta,'o')
grid on
axis([20 80 -1 1])
yline(-0.8);
I drew a line at y = -0.8, which is where I would suggest the baseline for your curve should roughly lie at.
That is, what is the asymptotic value of your function as x approaches infinity? You have written it as:
sin(theta) = -1+2*sqrt(alpha/x)*exp(-beta*(x-alpha)^2)
As x gets large, this approaches -1. Yet, clearly, your data does not get anywhere near -1, and it never will. Your data approaches a value closer to -0.8, MAYBE -0.85. So your data is inconsistent with your model. That you want to call the depedent variable sin(theta) does not matter. It simply will never get anywhere near -1.
Lets look at how you might fit this with the curve fitting toolbox. remember, that you will NEED good starting values for alpha and beta.
mdl = fittype('-1+2*sqrt(alpha./x).*exp(-beta*(x-alpha).^2)','indep','x');
fittedmdl = fit(x',sintheta',mdl,'startpoint',[20,0.01])
fittedmdl =
General model:
fittedmdl(x) = -1+2*sqrt(alpha./x).*exp(-beta*(x-alpha).^2)
Coefficients (with 95% confidence bounds):
alpha = 25.82 (23.66, 27.99)
beta = 0.01954 (0.0005392, 0.03854)
plot(fittedmdl)
hold on
plot(x,sintheta,'o')
yline(-.8);
Note that I estimated beta there, too, as the value for beta that you gave is also inconsistent with the data. If I try to force the model to use beta =1.135e-4;, I get garbage for a fit, far worse than what I got for the model I chose.
Now, look carefully at the curve that came out. It does not look at all like your data. The baseline is wrong. The shape of the curve in the peak is wrong.
Why do I say those shapes are "wrong"? The peak of the curve is fairly wide. It rolls over nice and slowly. And it never even manages to start coming back down, if I look at your data. But the transition region where it drops down from 1 towards -1, is VERY sharp. That means you need a relatively large value for beta, NOT one as small as you seem to think it should be. Otherwise, you cannot get a sharp transition. But with a sharp transition there, that also means the peak wants to drop down for small x, below the value of alpha. So your data is simply not consistent with that model.
I honestly don't know where you got that model. But it simply will never look like that data. Sorry. But trying to fit a square peg into a round hole tends only to damage either the peg or the hole.
(In fact, if I try to estimate the baseline from your data, it comes out to approximately -0.79, so my visual guess was pretty good.)

  0 Comments

Sign in to comment.


Answer by Matt J
on 9 Jun 2019
Edited by Matt J
on 10 Jun 2019

Are you sure the model isn't supposed to be a Gaussian+constant? It gives a much better fit. I used this FEX file to do the fit, but you could use the Curve Fitting toolbox as well, if you have it.
LB = [];
UB = {[],[],xdata(end)} ; %center the Gaussian left of the first data point
params0={-1,2,30,[]};
[params] = gaussfitn(xdata,ydata,params0,LB,UB);
[D,A,mu,sig]=deal(params{:})
fun=@(x) D + A*exp( -0.5 * (x-mu).' * inv(sig) *(x-mu) );
xc=linspace(xdata(1),xdata(end),1000);
plot(xdata,ydata,'o',xc,arrayfun(fun,xc),'-')

  0 Comments

Sign in to comment.