Python implementation of numerical calculation----Newton interpolation method

The biggest problem with Lagrange interpolation is that each time a new interpolation node is introduced, the base function changes, which is inappropriate in some actual production environments, and sometimes new measurement data is added to the interpolation node set.

Therefore, Newton interpolation method is formed by looking for the relationship between the interpolation function constructed by n interpolation nodes and the interpolation function constructed by n+1 interpolation nodes. The method of deducing Newton interpolation is an inductive method, that is, to compute Ln (x)-ln+1 (x), and to calculate the interpolation function of n+1 at the constant iteration from N=1.

The formula for Newton's interpolation method is:

Note: In the program I use W instead of

The key to calculate the Newton interpolation function is to calculate the difference quotient, the N-order difference quotient is represented as follows:

　　　　　　　　　　　　　　　　　　　　　　　　

I'm not talking about bad business here.

The formula for calculating the N-order difference quotient is this:

It is obvious that the N-order difference quotient needs to use two n-1, so it is easy to think of the recursive implementation to calculate the N-order difference quotient when programming, but it is necessary to pay attention to the potential danger of recursive stack overflow, when calculating the difference quotient

Even more so, each level of recursion contains two recursion, and the total number of recursion is full of two tree distribution:

　　　　

This means that the number of recursion will increase sharply: (. so in the specific application also need to change the idea or optimize the code according to the application .

Talk less, put the code over.

First, write the most critical step, that is, to calculate the N-Order difference quotient:

Look at the Newton interpolation function formula above, with the difference quotient, still poor

This is better achieved:

"" "@brief:  Get the WI (x) function;         WI means an example W1 = (x-x0); W2 = (x-x0) (X-X1); W3 = (x-x0) (x-x1) (x-x2) @param:  i-  i-order (i-polynomial) @param:  The horizontal axis of all the interpolated nodes @return: Returns the Wi (x) function "" "Def Get_wi (i = 0, Xi = []):    def Wi (x):        result = 1.0 for each in        range (i):            result *= (X-xi[each])        return result
   return Wi

　　　　

OK, Newton interpolation method The most important two parts have, the following is the two parts combination Cheng Nunton interpolation function, if it is C and other languages need to save some intermediate data, I used the Python closure directly return a Newton interpolation function, the closure can take advantage of its function in the context of the data.

"" "@brief: Get Newton interpolation function @" "" def get_newton_inter (xi = [], fi = []):    def newton_inter (x):        result = fi[0] for        i in Range (2, Len (xi)):            result + = (Get_order_diff_quot (xi[:i], fi[:i]) * GET_WI (I-1, xi) (x))        return result    Return Newton_inter

The above code is the translation of the Newton interpolation function formula, note that the parameters of the GET_WI function is i-1, the expression from the function can find the reason.

Construct some interpolation nodes

　

The interpolation node, where the interpolation node is generated with two functions, each of the two nodes x-axis distance of "    sr_x = [I for I in range ( -50, Wuyi, ten)]    sr_fx = [i**2 for i in sr_x]  

　　　

Get Newton interpolation function

    NX = Get_newton_inter (sr_x, Sr_fx)            # Get the interpolation function        tmp_x = [I for I in range ( -50, Wuyi)]          # test Case    tmp_y = [Nx (i) for I in tmp_x]               # Get the ordinate of the test case based on the interpolation function

　　

Full code:

#-*-Coding:utf-8-*-"" "Created on Thu Nov 18:34:21 2016@author:tete@brief: Newton interpolation Method" "" Import Matplotlib.pyplot as PLT "                                                        "@brief: Calculate N-order difference quotient f[x0, x1, x2 ... xn] @param: Xi all interpolation nodes ' horizontal axis collection                                                     O@param:fi the Ordinate collection/@return of all interpolated nodes: Returns the I-Order quotient of Xi (I is XI length minus 1)                                                        O o@notice:a. Must ensure that the XI is equal to the fi length                                           /\/b. Due to the use of recursion, so be careful not to explode the stack.                                                                                     o o o C. Recursive minus recursion (each layer recursion contains two recursive functions), each level of recursion is two power growth, the total number of times is a full two tree of all nodes (so very easy stack overflow) "" "Def get_order_diff_quot (xi = [], fi = []): If Len (xi) > 2 and Len ( FI) > 2:return (Get_order_diff_quot (Xi[:len (xi)-1], Fi[:len (FI)-1])-Get_order_diff_quot (Xi[1:len (xi)], fi [1:len (FI)])) /float (xi[0]-xi[-1]) return (fI[0]-fi[1])/float (xi[0]-xi[1]) "" @brief: Get the WI (x) function; WI means an example W1 = (x-x0); W2 = (x-x0) (X-X1); W3 = (x-x0) (x-x1) (x-x2) @param: I I-order (i-polynomial) @param: The horizontal axis of all the interpolated nodes @return: Returns the Wi (x) function "" "def get_wi (i = 0, Xi = [])    : Def Wi (x): result = 1.0 for all in range (i): Result *= (X-xi[each]) return result        return Wi "" "@brief: Get Newton interpolation function @" "" def get_newton_inter (xi = [], fi = []): Def newton_inter (x): result = Fi[0] for I in range (2, Len (xi)): Result + = (Get_order_diff_quot (xi[:i], fi[:i]) * GET_WI (I-1, xi) (x)) return result return Newton_inter "" "Demo:" "" "if __name__ = = ' __main__ ':" ' interpolation node, here two times The function generates an interpolation node, each two nodes x-axis distance of "sr_x = [I for I in range ( -50, Wuyi, ten)] Sr_fx = [i**2 for i in sr_x] Nx = GET_NEWT  On_inter (sr_x, SR_FX) # Get the interpolation function tmp_x = [I for I in range (-50, 51)] # test Case tmp_y = [Nx (i) for I in tmp_x] #According to the interpolation function to obtain the ordinate "' paint ' of the test case (" plt.figure ") Ax1 = Plt.subplot (111) Plt.sca (AX1) Plt.plot ( Sr_x, sr_fx, LineStyle = ', marker= ' o ', color= ' B ') Plt.plot (tmp_x, tmp_y, LineStyle = '--', color= ' R ') Plt.show ()

　　

　　

Python implementation of numerical calculation----Newton interpolation method