清空筛选条件
清空筛选条件

Fitting a hyperbola through a set of points

58 次查看（过去 30 天）
Ravi Mudragada
Ravi Mudragada 2022-9-14
评论： Ravi Mudragada 2022-9-15
I'm trying to fit a hyperbola through given data points. The code works fine for the first set of points but does not work for the other sets of points. The code based on this answer (https://in.mathworks.com/matlabcentral/answers/355943-how-can-we-fit-hyperbola-to-data#answer_280986) is below and the accompanying data is attached. Any suggestions would be highly encouraged!
clear all;
data=readtable("PI-No-StrainRate.xlsx");
x1=data.P1;
y1=data.I1;
%x1=data.P2; % Uncomment to replicate issue
%y1=data.I2; % Uncomment to replicate issue
hyprb = @(b,x1) b(1) + b(2)./(x1+b(3)); % Generalised Hyperbola
NRCF = @(b) norm(y1 - hyprb(b,x1)); % Residual Norm Cost Function
B0 = [0, 0, 0];
%options = optimset('MaxFunEvals',800000);
B1 = fminsearch(NRCF, B0);%options); % Estimate Parameters
figure(1)
%plot(y1, x1, '+r')
%x2 = [0:1:30];
hold on
plot(y1, x1, 'ko', hyprb(B1,x1), x1, '-r')
xlabel('Impulse (MPa-msec)')
ylabel('Pressure (MPa)')
hold on
%plot(y1,y1-hyprb(B1,x1))
grid on
text(15, 10, sprintf('y = %.4f %+.4f/(x %+.4f)', B1))
hold off
%x0 = [0,100];
%[B,resnorm,residual,exitflag,output] = lsqcurvefit(hyprb,B0,x1,y1);
%plot(x1,hyprb(B, x1),'b-')
%grid on

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2022-9-14
编辑：Torsten 2022-9-14
Ran in:
Remove the NaN value in the last position of data.P2 and data.I2:
clear all;
data=readtable("https://de.mathworks.com/matlabcentral/answers/uploaded_files/1124700/PI-No-StrainRate.xlsx");
%x1=data.P1;
%y1=data.I1;
x1=data.P2; % Uncomment to replicate issue
x1 = x1(1:end-1);
y1=data.I2; % Uncomment to replicate issue
y1 = y1(1:end-1);
hyprb = @(b,x1) b(1) + b(2)./(x1+b(3)); % Generalised Hyperbola
NRCF = @(b) norm(y1 - hyprb(b,x1)); % Residual Norm Cost Function
B0 = [0, 0, 0];
%options = optimset('MaxFunEvals',800000);
B1 = fminsearch(NRCF, B0);%options); % Estimate Parameters
figure(1)
%plot(y1, x1, '+r')
%x2 = [0:1:30];
hold on
plot(y1, x1, 'ko', hyprb(B1,x1), x1, '-r')
xlabel('Impulse (MPa-msec)')
ylabel('Pressure (MPa)')
hold on
%plot(y1,y1-hyprb(B1,x1))
grid on
text(15, 10, sprintf('y = %.4f %+.4f/(x %+.4f)', B1))
hold off
%x0 = [0,100];
%[B,resnorm,residual,exitflag,output] = lsqcurvefit(hyprb,B0,x1,y1);
%plot(x1,hyprb(B, x1),'b-')
%grid on

更多回答（1 个）

William Rose
William Rose 2022-9-14
Ran in:
@Ravi Mudragada,
The columns in the workbook are 12 long for P1, I1. The columns are 11 long for P2,I2.
When readtable() reads the workbook, it assigns NaN values to the final row of I2,P2. These Nan values are throwing off the fit. Do not include the last row of values in I2, P2, and it works fine.
data=readtable("PI-No-StrainRate.xlsx");
%x1=data.P1;
%y1=data.I1;
x1=data.P2(1:end-1); % Uncomment to replicate issue
y1=data.I2(1:end-1); % Uncomment to replicate issue
hyprb = @(b,x1) b(1) + b(2)./(x1+b(3)); % Generalised Hyperbola
NRCF = @(b) norm(y1 - hyprb(b,x1)); % Residual Norm Cost Function
B0 = [0, 0, 0];
%options = optimset('MaxFunEvals',800000);
B1 = fminsearch(NRCF, B0);%options); % Estimate Parameters
figure(1)
%plot(y1, x1, '+r')
%x2 = [0:1:30];
hold on
plot(y1, x1, 'ko', hyprb(B1,x1), x1, '-r')
xlabel('Impulse (MPa-msec)')
ylabel('Pressure (MPa)')
hold on
%plot(y1,y1-hyprb(B1,x1))
grid on
text(15, 10, sprintf('y = %.4f %+.4f/(x %+.4f)', B1))
hold off
Are you aware that you are plotting x on the vertical axis and y on the horizontal axis? That is up to you, of course.
I would define the fitting funciton as
hyprb = @(b,x1) b(1) + b(2)./(x1-b(3)); % note that I changed the sign on b(3)
because if you do, then b(1) is the coordinate of the vertical asymptote and b(3) is the coordinate of the horizontal asymptote, and with this change, the equation for the hyperbola can be rewritten
(y-b1)(x-b3)=b2
which has a nice symmetry. That version of the formula is not useful in your code, but it is nice for explaining the equation.
Good luck.
  3 个评论
显示 1更早的评论
William Rose
William Rose 2022-9-14
@Ravi Mudragada,
I did not see that @Torsten had posted his answer before I finished mine. I'm kind of slow and his was not up while I was composing mine. I recommend that you accept his since he delivered first.
Ravi Mudragada
Ravi Mudragada 2022-9-15
Did that! However your explanation of the issue helped a lot too. Much appreciated!

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 MATLAB 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by