Curve fitting with constraints
172 views (last 30 days)
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
Matt J on 1 Mar 2019
Edited: Matt J on 1 Mar 2019
Using the 'power1' fit type, I obtain pretty close agreement.
>> cfit = fit(x1(x1>0),y1(x1>0),'power1')
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)
More Answers (2)
John D'Errico on 1 Mar 2019
Edited: John D'Errico on 1 Mar 2019
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.
Alex Sha on 1 Feb 2020
Hellow, Paul, if adding one more parameter in your fitting function, i.e. y = a*exp(b*x) + c*exp(d*x)+f
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
Adjusted R-Square: 0.99999282087114
Determination Coef. (DC): 0.999991748082927
Parameter Initial Value