Curve fitting with constraints

Hi,

I am trying to fit an curve with an exponential function and need it to pass through (0,0) with a gradient of 0 at that point.

I've attached a picture of the data and curve. I am trying to use to the curve fitting toolbox to do this.

Any ideas would be greatly appreciated

Thanks

[Torsten]

Your conditions at 0 force the Exp2 function to be identically 0 - not a very interesting case for fitting:

f(0)=0 -> a+c = 0 -> c = -a

f'(0)=0 -> a*b + c*d = 0 -> a*b - a*d = 0 -> b = d

-> f(x) = a*exp(b*x) - a*exp(b*x) = 0

[Paul Van der Merwe]

Thanks,

I understand the premise, but I am struggling to apply to curve fitting in MATLAB.

[Torsten]

Why not choosing a parabola y=-a*x^n as fitting model ?

[Paul Van der Merwe]

A polynomial won't work. I get an error 'The equation is badly conditioned' and the fit for the exponential function is better. It just does not meet the two required conditions.

[Torsten]

Not a polynomial ; a and n are the paramters to be fitted in y=-a*x^n.

[Paul Van der Merwe]

The x values are negative and the fitted curve will look something like y=-a*x^1/n. I need the y values to be real, not complex.

[Torsten]

In your graphics, the x-values are positive (range from 0 to 90) and the y-values are negative( range from 0 to -3). If this is reversed in reality, you can fit the inverse function by the above ansatz.

[Paul Van der Merwe]

Sorry that's confusing. Ignore the axis labels (I need the inverted function). So the y values range of 0 - 90 and the x values range from 0 - -3e5. I am currently struggling to fit the inversed function with exponential terms.

Answer 1

[Matt J]

Using the 'power1' fit type, I obtain pretty close agreement.

>> cfit = fit(x1(x1>0),y1(x1>0),'power1')

cfit = 

     General model Power1:
     cfit(x) = a*x^b
     Coefficients (with 95% confidence bounds):
       a =    -0.07755  (-0.07888, -0.07622)
       b =       3.316  (3.312, 3.32)
       
>> plot(cfit,x1,y1)

[Paul Van der Merwe]

That's a great fit. However, I need the equation of that fit in terms of y (i.e. x = f(y)).

The inverted function from your reply will return complex numbers (as y has negative values) and I required real numbers. Any ideas?

[Matt J]

No, the inverted function

x=(y/a).^(1/n)

will not return complex numbers for negative y, because a is also negative.

[Paul Van der Merwe]

Very good point!

Thanks a lot for your help

[Nicholas Ross]

@Matt J thank you for posting this! I wasn't getting a good result with exp1 or lsqcurvefit and this fit my curve perfectly! Thanks for making my Monday morning

Answer 2

[Alex Sha]

Hellow, Paul, if adding one more parameter in your fitting function, i.e. y = a*exp(b*x) + c*exp(d*x)+f

Then:

from y(0) = 0, we get: f = -(a+c);

from y'(0) = 0, we get: d=-a*b/c;

It can be easy to get the:

Root of Mean Square Error (RMSE): 214.862401971985
Sum of Squared Residual: 434512996.964381
Correlation Coef. (R): 0.999996411191975
R-Square: 0.999992822396829
Adjusted R-Square: 0.99999282087114
Determination Coef. (DC): 0.999991748082927
Chi-Square: -26589.2854996605
F-Statistic: 569341000.214239

Parameter	Initial Value
a	-20088.9862898102
b	0.0320322385142336
c	55098.3960866006

Answer 3

[John D'Errico]

You cannot use the curve fitting toolbox to fit a model that has constraints on it like the ones you have. At least not directly.

Lets look at your model function however. The general model is:

y = a*exp(b*x) + c*exp(d*x)

Your constraints are that

y(0) = 0

y'(0) = 0

What do they mean in terms of the parameters?

y(0)=0 ==> a + b = 0

Or, we can write it as a = -b, essentially allowing us to eliminate b from the model completely. So now we can re-write the model in the form

y = a*exp(b*x) - a*exp(d*x)

Next, consider the slope boundary condition at x=0. Again, since exp(0) is 1, this reduces to:

y'(0) = 0 ==> a*b - a*c = 0

That is a problem however. Why? If we assume a is non-zero (as the problem becomes trivial if a is zero) then we now know that b == c. So we can reduce the model once more, into the form:

y(x) = a*exp(b*x) - a*exp(b*x) = 0

Does that seem problematic? It has a difference of two identical terms. So the ONLY exponential model that is a sum of exactly two exponentials of the form that you want, that also has y(0)=y-(0)=0, is the oddly trivial one:

y(x) = 0

So, can you solve this using the curve fitting toolbox? I suppose you might be able to do so, although I'm not sure how to enter a model that is constant at zero for all x.

I would strongly recommend you choose a model that does allow the required properties to hold, but also allows a fit to exist. Without seeing your data, it is very difficult to know what that model might be. For example, you might be able to use a very simple model of the form:

y(x) = a2*x.^2 + a3*x.^3 + a4*x.^4

etc. Don't go too high in the order of that model, as it will rapidly have numerical difficulties. You may need to consider scaling your data. A virtue of the model form I have suggested is it automatically forces y(0)=y'(0)=0.

My guess is that the data you have has relatively large x. So, how would I scale this problem? I would scale it as:

s = max(abs(x));

Now, solve the problem using the curve fitting toolbox by scaling x and y as:

xs = x/s;
ys = y/s^2;
ft = fittype('a2*x.^2 + a3*x.^3 + a4*x.^4','indep','x');
mdl = fit(xs(:),ys(:),ft);

Don't forget to unscale the coefficients afterwards. If you need a better recommendation for a model, you need to provide the data.

[Paul Van der Merwe]

Thank you, that makes a lot of sense and explains it well. However, using a polynomial as you suggest causes a warning 'Equation is badly conditioned.'

I have attached the data points, if that helps. I am trying to fit a curve to the inverted function (i.e. x = f(y)) with the boundary conditions.

An exponential fit seems the best method as I need real numbers, not complex.

Thanks for explaining so far!