How to adapt a function to a change of two points?

Say you have a function defined in Excel, and you fit the data to a polynomial curve using cftool. Now you want to change just two points of the function, from (x,y)_original to (x,y)_new. You want the rest of the function to adapt to that change. How do you do this?

Something like this:

Right now I have a rational fit with 9 coefficients (5th power numerator and 5th power denominator, 5 coefficients in the numerator and 4 coefficients in the denominator while the coefficient of x^5 in the denominator is 1). I can, however, choose a different fit (polynomial, sum of sines etc.) if it would help here.

I have tried changing two coefficients at random to force the function to fit the two new points, but it produces something that doesn't make sense.

Here's what I have:

This is a rational fit to a set of Excel data produced with Matlab cftool:

```

p1 = 0.09503;

p2 = 0.8033;

p3 = 2.654;

p4 = 4.375;

p5 = 3.691;

p6 = 1.221;

q1 = 6.399;

q2 = 21.65;

q3 = 34.98;

q4 = 29.02;

q5 = 9.704;

m = (x - 0.9857)./0.1235;

y = (p1*m.^5 + p2*m.^4 + p3*m.^3 + p4*m.^2 + p5*m + p6) ./ (m.^5 + q1*m.^4 + q2*m.^3 + q3*m.^2 + q4*m + q5);

```

I need to make the function go through points (0.78; 0.002) and (1.05473536553; 9.869964329717e-2), while keeping the shape.

Answer 1

[Mathieu NOE]

hello Daniel

this is my suggestion - a simple anamorphosis of the curve along the vertical axis (= a vertical shift and a multiplicative correction)

plot :

code

x = 0.6:0.01:1.2;

p1 = 0.09503;
p2 = 0.8033;
p3 = 2.654;
p4 = 4.375;
p5 = 3.691; 
p6 = 1.221;  
q1 = 6.399;
q2 = 21.65;  
q3 = 34.98;  
q4 = 29.02; 
q5 = 9.704;  
m = (x - 0.9857)./0.1235;
y = (p1*m.^5 + p2*m.^4 + p3*m.^3 + p4*m.^2 + p5*m + p6) ./ (m.^5 + q1*m.^4 + q2*m.^3 + q3*m.^2 + q4*m + q5);


% I need to make the function go through points A = (0.78; 0.002) 
% and B  = (1.05473536553; 9.869964329717e-2), while keeping the shape.

% what are the y coordinates for x = 0.78 and x = 1.05473536553 ? 
x1 = 0.78;
m = (x1 - 0.9857)./0.1235;
y1 = (p1*m.^5 + p2*m.^4 + p3*m.^3 + p4*m.^2 + p5*m + p6) ./ (m.^5 + q1*m.^4 + q2*m.^3 + q3*m.^2 + q4*m + q5); % current  coordinate
y1b = 0.002; % new target coordinate

x2 = 1.05473536553;
m = (x2 - 0.9857)./0.1235;
y2 = (p1*m.^5 + p2*m.^4 + p3*m.^3 + p4*m.^2 + p5*m + p6) ./ (m.^5 + q1*m.^4 + q2*m.^3 + q3*m.^2 + q4*m + q5); % current  coordinate
y2b = 9.869964329717e-2; % new target coordinate


% solution : search a, b coeff so that : 
% y1b = a * y1 + b
% y2b = a * y2 + b

a = (y1b - y2b) / (y1 - y2);
b =  y1b - a * y1;

% plot new line 
ynew = a * y + b;

plot(x,y,'b',x1,y1b,'*r',x2,y2b,'*k',x,ynew,'g.-');
legend(' original line','new point A','new point B',' new line');

[Daniel Powell]

Hi Mathieu. Great answer! I forgot to add that the second point is meant to be a maximum. How would I implement this into your code? Also, how could I get rid of the imperfection mid-curve? Thanks.

[Mathieu NOE]

Hi Daniel

the problem with the mid curve imperfection can be shown with greater accuracy if I refine the x axis vector ; with a delta x = 1e-3 , the mid curve imperfection looks like this

this is the cause of having a denominator in your model which can creates poles (the peak / dip is due to that) - why not stick whith for example a simple polynomial model (no denominator)

I tried with an order 16 , and to not have the imperfection I generated the data with a coarse x resolution (trying not to have a x value close to the "danger zone"), then refined the output with spline interpolation. We can see that a 16th order polynomial can fit quite well the shape of the curve

so this is just another topic, but it may as well be of some interest for you :

clc
clearvars

x = 0.6:0.03:1.2;       % to hide the imperfection mid-curve
% x = 0.6:0.001:1.2;       % to show the imperfection mid-curve

p1 = 0.09503;
p2 = 0.8033;
p3 = 2.654;
p4 = 4.375;
p5 = 3.691; 
p6 = 1.221;  
q1 = 6.399;
q2 = 21.65;  
q3 = 34.98;  
q4 = 29.02; 
q5 = 9.704;  
m = (x - 0.9857)./0.1235;
y = (p1*m.^5 + p2*m.^4 + p3*m.^3 + p4*m.^2 + p5*m + p6) ./ (m.^5 + q1*m.^4 + q2*m.^3 + q3*m.^2 + q4*m + q5);
% I need to make the function go through points A = (0.78; 0.002) 
% and B  = (1.05473536553; 9.869964329717e-2), while keeping the shape.

% spline interpolation to get higher resolution plot (without the
% imperfection)
xx = linspace(min(x),max(x),200);
yy = interp1(x,y,xx,'spline');

% figure(1);
% plot(x,y,'-',xx,yy,'-')
% legend('raw data','interpolated')

% simple polyfit without denominator
% Fit a polynomial p of degree "degree" to the (x,y) data:
degree = 16;
p = polyfit(xx,yy,degree);
% Evaluate the fitted polynomial p and plot:
f = polyval(p,xx);

figure(2);plot(x,y,'*b',xx,yy,'-',xx,f,'--')
legend('raw data','interpolated','polyfit')

[Mathieu NOE]

hello again

so this is the full code with vertical + horizontal corrections as your requested

clc
clearvars

x = 0.6:0.03:1.2;       % to hide the imperfection mid-curve
% x = 0.6:0.001:1.2;       % to show the imperfection mid-curve

p1 = 0.09503;
p2 = 0.8033;
p3 = 2.654;
p4 = 4.375;
p5 = 3.691; 
p6 = 1.221;  
q1 = 6.399;
q2 = 21.65;  
q3 = 34.98;  
q4 = 29.02; 
q5 = 9.704;  
m = (x - 0.9857)./0.1235;
y = (p1*m.^5 + p2*m.^4 + p3*m.^3 + p4*m.^2 + p5*m + p6) ./ (m.^5 + q1*m.^4 + q2*m.^3 + q3*m.^2 + q4*m + q5);
% I need to make the function go through points A = (0.78; 0.002) 
% and B  = (1.05473536553; 9.869964329717e-2), while keeping the shape.

% spline interpolation to get higher resolution plot (without the
% imperfection)
xx = linspace(min(x),max(x),200);
yy = interp1(x,y,xx,'spline');

% simple polyfit (without denominator)
% Fit a polynomial p of degree "degree" to the (x,y) data:
degree = 16;
p = polyfit(xx,yy,degree);
% Evaluate the fitted polynomial p and plot:
f = polyval(p,xx);

figure(1);plot(x,y,'*b',xx,yy,'-',xx,f,'--')
legend('raw data','interpolated','polyfit')

%% step 1 : vertical axis correction 
% what are the y coordinates for x = 0.78 and x = 1.05473536553 ? 
x1 = 0.78;
m = (x1 - 0.9857)./0.1235;
y1 = (p1*m.^5 + p2*m.^4 + p3*m.^3 + p4*m.^2 + p5*m + p6) ./ (m.^5 + q1*m.^4 + q2*m.^3 + q3*m.^2 + q4*m + q5); % current  coordinate
y1b = 0.002; % new target coordinate
x2 = 1.05473536553;
m = (x2 - 0.9857)./0.1235;
y2 = (p1*m.^5 + p2*m.^4 + p3*m.^3 + p4*m.^2 + p5*m + p6) ./ (m.^5 + q1*m.^4 + q2*m.^3 + q3*m.^2 + q4*m + q5); % current  coordinate
y2b = 9.869964329717e-2; % new target coordinate

% solution : search a, b coeff so that : 
% y1b = a * y1 + b
% y2b = a * y2 + b
a = (y1b - y2b) / (y1 - y2);
b =  y1b - a * y1;

% plot new line 
ynew = a * yy + b;
% figure(2)
% plot(xx,yy,'b',x1,y1b,'*r',x2,y2b,'*k',xx,ynew,'g.-');
% legend(' original line','new point A','new point B',' new line');

%% step 2 : horizontal axis correction 
[val, ind] = max(ynew);
x_max = xx(ind);

% solution : search a, b coeff so that : 
% x1 = a * x1 + b
% x2 = a * x_max + b
a = (x1 - x2) / (x1 - x_max);
b =  x1 - a * x1;

% plot new line 
xnew = a * xx + b;
ind = find(xnew>=min(x) & xnew<=max(x));

figure(2)
plot(xx,yy,'b',x1,y1b,'*r',x2,y2b,'*k',xx,ynew,'g--',xnew(ind),ynew(ind),'k.-');
legend(' original line','new point A','new point B',' new line 1 (vert. correction)',' new line 2 (vert. + hor. correction)');

[Daniel Powell]

Hi Mathieu. Fantastic answer! You've completely solved the question I had, thank you.

[Mathieu NOE]

My pleasure !

Would you mind accepting my answer ?

tx