1 view (last 30 days)

Show older comments

I currently have a code where I am initializing many parameters, that will be solved for later.

I do this by:

C1=params(1);

C2=params(2);

C3=params(3);

C4=params(4);

C5=params(5);

C6=params(6);

C7=params(7);

C8=params(8);

C9=params(9);

C10=params(10);

C11=params(11);

C12=params(12);

C13=params(13);

There's got to be an easier way...

Additionally, I set my intial guesses like:

params0=[0.1,0.1,0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1];

Is there a way to more easily set this, assuming they all have the same initial guess? Typing them all out that way is tedious, and my code is constantly changing. How can I automate it?

Matt J
on 19 Mar 2019

Edited: Matt J
on 19 Mar 2019

Is there a way to more easily set this, assuming they all have the same initial guess?

One way:

params0=linspace(0.1,0.1,13)

Another is:

params0=repelem(0.1,1,13);

There's got to be an easier way...

If you have a very large vector of variables, it does not make sense to assign them to separate variables. You should probably be manipulating them as a vector. But if you must index the variables individually, I would just rename params to something shorter,

C=params;

and then refer to the variables via indexing C(1), C(2),..C(13) throughout the code. Although, it is then unclear why you chose to name the vector params in the first place...

Matt J
on 20 Mar 2019

Edited: Matt J
on 20 Mar 2019

One other thing I would point out is the model function you have shown us here is linear in the coefficients. This means that it can be represented as a matrix multiplication and the entire fit can be done with simple linear algebra instead of an iterative solver like lsqcurvefit.

A=func2mat(@(p) modelfun(p,xdata), params0);

coefficients =A\g(:);

That's it!

Matt J
on 20 Mar 2019

Just to sum things up here, below is my implementation of the fit, combining all my various recommendations. Both the algebraic method that I proposed and the method using lsqcurvefit are implemented and compared. You will see that the algebraic method is both faster and more accurate (gives a lower resnorm) than lsqcurvefit.

%% Data Set-up

num=xlsread('example.xlsx',1,'A2:F18110');

Np=4; %polynomial order

g = num(:,6);

otherStuff.r = num(:,4);

otherStuff.eta = num(:,3);

otherStuff.g = g;

otherStuff.eta_c = 0.6452;

otherStuff.Np=Np;

otherStuff.A1 = -0.4194;

otherStuff.A2 = 0.5812;

otherStuff.A3 = 0.6439;

otherStuff.A4 = 0.4730;

%% Fitting (2 methods)

%fit using matrix algebra

tic;

A=func2mat(@(p) modelfun(p,otherStuff), ones(Np+1,Np+1));

Cfit1 =A\g(:);

resnorm1=norm(A*Cfit1-g(:));

toc %Elapsed time is 0.124593 seconds.

%fit iteratively with lsqcurvefit

tic;

options=optimoptions(@lsqcurvefit,'MaxFunctionEvaluations',10000);

Cinitial(1:Np+1,1:Np+1)=0.1;

[Cfit2, resnorm2]=lsqcurvefit(@modelfun,Cinitial,otherStuff,g,[],[],options);

toc %Elapsed time is 1.687990 seconds.

%% Display

figure(1); showFit(Cfit1,otherStuff);

figure(2); showFit(Cfit2,otherStuff);

resnorm1,resnorm2,

function ghat=modelfun(C,otherStuff)

r = otherStuff.r;

eta = otherStuff.eta;

eta_c = otherStuff.eta_c;

Np = otherStuff.Np;

A1 = otherStuff.A1;

A2 = otherStuff.A2;

A3 = otherStuff.A3;

A4 = otherStuff.A4;

PV = 1 + 3*eta./(eta_c-eta)+ A1*(eta./eta_c) + 2*A2*(eta./eta_c).^2 +...

3*A3*(eta./eta_c).^3 + 4*A4*(eta./eta_c).^4;

rdf_contact = (PV - 1)./(4*eta);

Cflip=rot90(reshape(C,Np+1,Np+1),2);

poly_guess = polyVal2D(Cflip,r-1,eta/eta_c,Np,Np);

ghat = (poly_guess.*rdf_contact);

end

function showFit(Cfit,otherStuff)

r = otherStuff.r(:);

eta = otherStuff.eta(:);

g = otherStuff.g(:);

ghat=modelfun(Cfit,otherStuff);

tri = delaunay(r,eta);

%% Plot it with TRISURF

h=trisurf(tri, r, eta, ghat);

h.EdgeColor='b';

h.FaceColor='b';

axis vis3d

hold on;

scatter3(r,eta,g,'MarkerEdgeColor','none',...

'MarkerFaceColor','r','MarkerFaceAlpha',.05);

xlabel 'r', ylabel '\eta', zlabel 'g(r,\eta)'

hold off

end

Both methods give a similar looking surface plot:

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!