Curve fitting with custom function
显示 更早的评论
Hello community,
I am trying to do a curve fitting on some experimental data with a custom function. I am more specifically trying to bring closer my function to the data. The function is the following one:
Every parameters of the function, except the shear rate (
) , can vary in order to fit best the data.
) , can vary in order to fit best the data. Here is a graph with the data and the function plotted with initial values:
The initial values are the following ones:
muInf=0.0035 %[Pa.s]
mu0=0.108; %[Pa.s]
lambda=8.2;
n=0.3;
a=0.64;
viscCar=muInf+(mu0-muInf)*(1+(lambda.*shearRate).^a).^((n-1)/a);
loglog(shearRate,viscCar,'k');
I already tried to use the curve fitting toolbox but i wasn't able to keep the look of the function.
Does anyone know how i can do that ?
Thank you for the help !
采纳的回答
It would help to have the data.
shearRate = logspace(-3, 4, 100); % Create Data
muv = 0.5-tanh(shearRate)*0.01 + randn(size(shearRate))*0.01; % Create Data
shearRate = shearRate(:);
muv = muv(:);
% muInf=0.0035; %[Pa.s]
% mu0=0.108; %[Pa.s]
% lambda=8.2;
% n=0.3;
% a=0.64;
B0 = [0.0035; 0.108; 8.2; 0.3; 0.64];
viscCar = @(muInf,mu0,lambda,a,n,shearRate) muInf+(mu0-muInf)./(1+(lambda.*shearRate).^a).^((n-1)/a);
viscCarfcn = @(b,shearRate) viscCar(b(1),b(2),b(3),b(4),b(5),shearRate);
mumdl = fitnlm(shearRate,muv, viscCarfcn, B0)
Beta = mumdl.Coefficients.Estimate
figure
loglog(shearRate,viscCarfcn(Beta,shearRate),'k');
hold on
plot(shearRate, muv, '.b')
hold off
grid
This works, and would actually make sense with the actual data.
.
4 个评论
Da125
2021-11-17
Thank you Star Strider for your answer !
Yes i attach the .csv file with the data.
I tried to implement your suggestion but i have issues with that. I obtained some negative coefficient that don't make sense. Here is what i obtained:
With the coefficients:
There were several errors in ‘viscCar’ that I found and coirrected. The most important was replacing the vectorised multiplication of the two terms by a vectorised division, and the others were adding appropriate parentheses to the denominator and replacing ‘(n-1)’ with ‘(1-n)’ in the expoenet. With those corrections, it works. (I added the fminsearch call to estimate ‘B0’ for fitnlm because fminsearch uses a derivative-free search, and would likely find a suitable set of parameters relatively quicklly that it can then provide to fitnlm.)
The principal problem, however, was mis-coding the ‘visccar’ objective function.
T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/803909/bloodExpData.csv', 'VariableNamingRule','preserve');
shearRate = T1.Var1;
muv = T1.Var2;
% shearRate = logspace(-3, 4, 100); % Create Data
% muv = 0.5-tanh(shearRate)*0.01 + randn(size(shearRate))*0.01; % Create Data
shearRate = shearRate(:);
muv = muv(:);
% muInf=0.0035; %[Pa.s]
% mu0=0.108; %[Pa.s]
% lambda=8.2;
% n=0.3;
% a=0.64;
% B0 = [0.0035; 0.108; 8.2; 0.3; 0.64];
B0 = rand(1,5);
viscCar = @(muInf,mu0,lambda,a,n,shearRate) muInf+(mu0-muInf)./((1+(lambda.*shearRate).^a).^((1-n)/a));
viscCarfcn = @(b,shearRate) viscCar(b(1),b(2),b(3),b(4),b(5),shearRate);
[B0,fv] = fminsearch(@(b) norm(viscCarfcn(b,shearRate) - muv), B0)
mumdl = fitnlm(shearRate,muv, viscCarfcn, B0)
Beta = mumdl.Coefficients.Estimate
figure
loglog(shearRate,viscCarfcn(Beta,shearRate),'k');
hold on
plot(shearRate, muv, '.b')
hold off
grid
Success!
.
Da125
2021-11-18
Star Strider
2021-11-18
As always, my pleasure!
.
更多回答（1 个）
Alex Sha
2021-11-18
It is hard to get stable and unique result for Da125's problem, especially for parameters of "lambda" and "n". refer to the result below:
Root of Mean Square Error (RMSE): 0.00032178294087402
Sum of Squared Residual: 2.89923930905092E-6
Correlation Coef. (R): 0.998376698597668
R-Square: 0.996756032302778
Parameter Best Estimate
---------- -------------
muinf -0.0902678665303079
mu0 0.00320052355322768
lambda 207505339.640542
a -0.504386079091155
n -1252.71593136491
1 个评论
Star Strider
2021-11-18
.
类别
在 帮助中心 和 File Exchange 中查找有关 Get Started with Curve Fitting Toolbox 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
