9 views (last 30 days)

Show older comments

I have a code that is of a function g(X,Y). The below code fits my data very well. However, I want to add some constraints at the boundaries where my data is not. For instance, how can I set the coefficient such that when Y=0, then the value of my function is 1?

Can anyone help with this?

clear

clc

close all

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

eta = num(:,3);

r = num(:,4);

g = num(:,6);

%Do the surface fit

options=optimoptions(@lsqcurvefit,'SpecifyObjectiveGradient',false,'MaxFunctionEvaluations',50000,'MaxIterations',1000);

xdata={r,eta};

params0=linspace(0.1, 0.1, 50);

[params, resnorm]=lsqcurvefit(@modelfun,params0,xdata,g,[],[],options)

for i=1:50

C(i)=params(i);

end

xmin = 1; xmax = 2; dx = 0.01;

etamin=0; etamax=0.55; deta=0.01;

[Xgrid,etagrid]=ndgrid(xmin:dx:xmax,etamin:deta:etamax);

surf(Xgrid,etagrid,modelfun(params,{Xgrid,etagrid}))

zlim([0 8]);

hold on;

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

'MarkerFaceColor','b','MarkerFaceAlpha',.5); hold off

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

function [out,Jacobian]=modelfun(params,xdata)

X=xdata{1};

Y=xdata{2};

for i=1:50

C(i)=params(i);

end

A1 = -0.4194;

A2 = 0.5812;

A3 = 0.6439;

A4 = 0.4730;

eta_c = 0.6452;

PV = 1 + 3*Y./(eta_c-Y)+ A1*(Y./eta_c) + 2*A2*(Y./eta_c).^2 + 3*A3*(Y./eta_c).^3 + 4*A4*(Y./eta_c).^4;

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

poly_guess = polyVal2D(C,X-1,Y/eta_c,4,4);

out = (poly_guess.*rdf_contact);

if nargout>1

%Jacobian=[X(:), Y(:), ones(size(X(:)))];

end

end

Adam Danz
on 20 Mar 2019

"...I want to add some constraints at the boundaries... "

In your call to lsqcurvefit(), you're not using the upper and lower bounds options (5th and 6th inputs are empty). I'd start out by imposing boundaries.

"For instance, how can I set the coefficient such that when Y=0, then the value of my function is 1"

I don't follow this part. The variable "Y" doesn't exist in your sample code. Also, what do you mean "the value of my function"? Do you mean the output of your nonlinear function or do you mean the coefficient estimates produced by the fit?

Matt J
on 21 Mar 2019

Edited: Matt J
on 21 Mar 2019

Here is an adaptation of the algebraic solution that I presented in your previous thread. I don't bother with the lsqcurvefit alternative, since as I showed you, it is both slower and less accurate. Below are also surface plots of the fitted surface, both in the neighborhood of the data points and near eta=0 where you can see that everything converges to the target value parameter, g0, as eta-->0. In this example, I have set g0=1.3.

clear

clc

close all

%% Data Set-up

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

Np=4; %polynomial order

g0=1.3; %target value at eta=0

g = num(:,6);

r = num(:,4);

eta = num(:,3);

otherStuff.r = r;

otherStuff.eta = eta;

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;

%% Fit using matrix algebra

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

q=4*g0./((3+otherStuff.A1)/otherStuff.eta_c); %weight on limiting target value as eta-->0

b=g(:)-q*A0(:,1);

A=A0;

A(:,1:Np+1)=[];

coeffs =[zeros(Np+1,1); A\b]; %fitted coefficients

coeffs(1)=q;

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

%% Also check surface near eta=0

[Rg,Eg]=ndgrid(unique(r),linspace(1e-8,0.01,100));

differentStuff=otherStuff;

differentStuff.r = Rg(:);

differentStuff.eta = Eg(:);

differentStuff.g=nan([numel(Rg),1]);

figure(2); showFit(coeffs,differentStuff);

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(coeffs,otherStuff)

r = otherStuff.r(:);

eta = otherStuff.eta(:);

g = otherStuff.g(:);

ghat=modelfun(coeffs,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

Matt J
on 21 Mar 2019

Any chance you could modify the code without the weighting function and with the new data? Note that this means rdf_contact and all the coeffs (A1, A2, ..) with it can be removed.

It's virtually trivial now, with the Curve Fitting Toolbox:

%% Load Data

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

eta_c= 0.6452;

r = num(:,4);

eta = num(:,3);

H = num(:,5);

%% Set-up for fit

[I,J]=ndgrid(0:4);

Terms= compose('C%u%u*r^%u*eta^%u', I(:),J(:),I(:),J(:)) ;

Terms=strjoin(Terms,' + ');

independent={'r','eta'};

dependent='H';

knownCoeffs= {'C00','C10','C20','C30','C40', 'C01','C11','C21','C31','C41'};

knownVals=num2cell([ [1 , 0 , 0 , 0 , 0 ], [ 8 , -6 , 0 , 0.5 , 0 ]/eta_c ]);

allCoeffs=compose('C%u%u',I(:),J(:));

[unknownCoeffs,include]=setdiff(allCoeffs,knownCoeffs);

ft=fittype(Terms,'independent',independent, 'dependent', dependent,...

'coefficients', unknownCoeffs,'problem', knownCoeffs);

%% Fit and display

fobj = fit([r,eta],H,ft,'problem',knownVals) ;

hp=plot(fobj,[r,eta],H);

set(hp(1),'FaceAlpha',0.5);

set(hp(2),'MarkerFaceColor','r');

xlabel 'r',

ylabel '\eta'

zlabel 'H(r,\eta)'

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

Start Hunting!