Polynomial fitting going through 1 point with force derivative =0

Question

I am looking for something similar to numpy.polyfit I have a curve with a fixed point that looks like a degree-2 polynomial

What i want to do is :

going through the first point exactly (in the following example (0.05 , 1.0))
On the first point to have a derivative =0

example :

TabX :

0,050 ; 0,055 ; 0,060 ; 0,065 ; 0,070 ; 0,075 ; 0,080 ; 0,085 ; 0,090 ; 0,095 ;     0,100 ; 0,110 ; 0,120 ; 0,130 ; 0,140 ; 0,150 ; 0,160 ; 0,170 ; 0,180 ; 0,190 ; 0,200 ; 0,210 ; 0,220 ; 0,230 ; 0,243

TabY :

1,000000000 ; 1,000446069 ; 1,000395689 ; 1,000359466 ; 1,000313867 ; 0,999937484 ; 0,999760969 ; 0,999811533 ; 0,999966352 ; 0,999767956 ; 1,000148567 ; 1,000634904 ; 1,000735849 ; 1,001199937 ; 1,001510678 ; 1,001722496 ; 1,001992602 ; 1,002487029 ; 1,003492247 ; 1,004006533 ; 1,004832845 ; 1,005730132 ; 1,006327527 ; 1,007109894 ; 1,008266254

I already found a "simple but barbaric" solution for going through the first point : I add a lot of weight to this point, either with weight function, or by adding a lot of 0.05 in TabX and a lot of 1.0 in TabY and using the normal np.polyfit function. It is ugly but it works.

But i really do not know how to had the derivative=0 in the point (0.05;1.0)

also, in this thread it is said that Lagrange multiplier could do the trick, but i was not able to use the function written by @Jaime, i have the raised error LinAlgError, 'Singular matrix'

Also, as I have to be able to launch this script in a basic Abaqus, the solution can only use basic python 2.7 and numpy. There is no way to use scipy ormatplotlib within this implementation.

_____________________________________________________________________

Edit : Using DanielF answer, I could make something halfworking (also thank you DanielF for your correction on my original post, much simpler to read)

Using a new origin is a good idea however, i cannot make it efficient with my data this is working :

def WorkingDeg4 ():
    x = np.arange(100)
    y0 = 4.0*x**4+0.07 * x ** 3 + 0.3 * x ** 2 + 1.1 * x
    y = y0 + 1000 * np.random.randn(x.shape[0])
    XX = np.vstack((x**4,x ** 3, x ** 2, x, np.ones_like(x))).T
    p_all = np.linalg.lstsq(XX, y)[0]
    pp = np.polyfit(x, y, 3)
    p_no_offset = np.linalg.lstsq(XX[:, :-1], y)[0]
    y_fit = np.dot(p_no_offset, XX[:, :-1].T)
    for i in range(0,len(x)):
        print x[i],y0[i],y[i],y_fit[i]

But if i want to make int using my data in the main i put :

if __name__ == '__main__':
    MyDataX=[0.050 , 0.055 , 0.060 , [...], 0.243]
    MyDataY=[1.000000000 , 1.000446069 , 1.000395689 , [...] , 1.008266254]
    TabX=[0.0]*len(MyDataX)
    TabY=[0.0]*len(MyDataY)
    for i in range(0,len(TabX)):
        TabX[i]=MyDataX[i]-MyDataX[0]
        TabY[i]=MyDataY[i]-MyDataY[0]

So, at this point, I did the "back to origin" phase

And i want to do the same as the def Working but for my data, so I did a copy past of tge WorkingDeg4 and just go rid of the creation of x and y and put it into argument

def NOTWorkingDeg4 (x,y): 
    XX = np.vstack((x**4,x ** 3, x ** 2, x, np.ones_like(x))).T
    p_all = np.linalg.lstsq(XX, y)[0]
    pp = np.polyfit(x, y, 3)
    p_no_offset = np.linalg.lstsq(XX[:, :-1], y)[0]
    y_fit = np.dot(p_no_offset, XX[:, :-1].T)

and this one is not working .... I have a mystake at the line

XX = np.vstack((x**4,x ** 3, x ** 2, x, np.ones_like(x))).T

saying TypeError: unsupported operand type(s) for ** or pow(): 'list' and 'int'

So for what i understand, it does not want to do it because my when it does x**4 , x is not an integer. But i do not know hot to go around the problem

_____________________________________________________________________ edit 2 : found it : the problem was to initiate TabX and TabY, not as array, but as np.array wrong :

TabX=[0.0]*len(MyDataX)
TabY=[0.0]*len(MyDataY)

Right :

TabX=np.array([0.0]*len(LongueursFissureGlobale))
TabY=np.array([0.0]*len(CourbeInterpolationGlobale))

You are probably going to want to shift your data so that `(0.5, 1.0)` is your origin `(0, 0)` to make your calculations easier, and then do a `numpy.linalg.lstsq`. See [This answer](https://stackoverflow.com/questions/32260204/numpy-polyfit-passing-through-0). You can then set the first derivative to zero by not only dropping the constants from the `a` matrix, but the `x` values as well. — Daniel F, Aug 31 '17 at 13:27
if you look at how a minimization of this sort is computed, it should be pretty easy to derive the modified formula. Will write something up in a few hours. — Mad Physicist, Aug 31 '17 at 20:23

score 1 · Accepted Answer · answered Aug 31 '17 at 13:41

1

You are probably going to want to shift your data so that (0.05, 1.0) is your origin (0, 0) to make your calculations easier, and then do a numpy.linalg.lstsq

x = TabX - 0.05
y = TabY - 1.0
X_poly = np.vstack((x ** 4, x ** 3, x ** 2))
poly_coeffs = np.linalg.lstsq(X_poly.T, y)
y_fit = np.dot(poly_coeffs, X_poly)

If it's needed, you'll have to transform the polynomial back to your old coordinates, but this transform makes the fitting much simpler.

See This answer for more details.

answered Aug 31 '17 at 13:41

Daniel F

13,620
2
29
55

Yeah, when dealing with polynomial derivatives, setting as many things to zero as possible makes calculations much easier. If this answers your question remember to mark the check. – Daniel F Sep 01 '17 at 12:54
It was one step, but i needed time to edit the first post since i could not resolve it yet. But thank you for your help. – texas Sep 01 '17 at 13:08

Polynomial fitting going through 1 point with force derivative =0

1 Answers1