6 views (last 30 days)

Show older comments

Hello,

I have the following code, which produces a surface of my data. (Example xlsx file attached). Certain coeffs are already constrained. My question is, how do I constrain the fit so that H(r,eta) does not go above 1? Also, even though I do not have data past eta/eta_c of around 0.15, how do I extend my fitted surface to eta/eta_c = 0? At 0, the function should equal 1 since the C00 coefficient is set to 1.

close all

%% Load Data

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

eta_c= 0.6452;

r = num(:,4);

eta = num(:,3);

H = num(:,8);

%% Set-up for fit

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

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/eta_c],H,ft,'problem',knownVals);

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

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

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

xlabel 'r',

ylabel '\eta/\eta_c'

zlabel 'H(r,\eta)'

zlim([0 1.1]);

ylim([0 0.55]);

Matt J
on 25 Mar 2019

Edited: Matt J
on 25 Mar 2019

My question is, how do I constrain the fit so that H(r,eta) does not go above 1?

Polynomials cannot be bounded everywhere. You would at the very least have to decide on a particular region where this bound would apply.

At 0, the function should equal 1 since the C00 coefficient is set to 1.

For that, you need to fix this line:

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

how do I extend my fitted surface to eta/eta_c = 0

Just pass extra points to the plot command,

extraR=[min(r);max(r)]; extraEta=[0;0]; extraH=fobj(extraR,extraEta/eta_c);

Rp=[extraR;r]; Etap=[extraEta; eta]/eta_c; Hp=[extraH;H] ;

hp=plot(fobj,[Rp,Etap],Hp);

Matt J
on 28 Mar 2019

Hmmm. Simpler than I thought.

Etas=0.2:0.2:0.8 ;

for i=1:numel(Etas)

f_restricted = @(r) fobj(r(:).', Etas(i)/eta_c ) ;

figure(i), fplot(f_restricted)

end

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

Start Hunting!