Implementation and analysis of polynomial interpolation algorithm based on numerical calculation algorithm

Numerical calculation refers to the algorithm in the field of numerical analysis. Numerical analysis is a special study of data problems related to numerical and approximate values, and numerical calculation plays a particularly important role in the study of numeric analysis.

polynomial interpolation is a method of calculating the approximate value of a function . Where the function value is known only on a few points.

The basis of this algorithm is to establish an interpolation polynomial pn (z) with a progression of less than or equal to N, where n+1 is the number of points known as the value of the function.

polynomial interpolation method

Many problems can be described in terms of functions. However, this function is usually unknown, and we can infer the approximate model of the function by a small number of known points. For this purpose, interpolation is done between known points. , about function f (x), the known point is x0...x8, represented by a black dot in the graph. With the help of interpolation, we are able to get the values of functions at Z0, Z1, Z2, and in small white squares in the graph. This section mainly discusses the polynomial interpolation method.

the fundamental point of polynomial interpolation is to establish a special form of polynomial, called interpolation polynomial .

To understand the meaning of interpolation polynomial in depth, let's look at some basic laws of the polynomial:

First, the polynomial is a function with the following form:

P (x) = a0 + a1x + a2x2 + ... + anxn

The A0,...,an here is the coefficient . when an is a nonzero integer, this form of polynomial is called an N-order polynomial . This is the exponential form of the polynomial, which is especially common in mathematical problems. However, other forms of polynomial are easier in some specific environments. For example, in the polynomial interpolation problem, Newton interpolation polynomial is a good example :

P (x) = a0 + A1 (X-C1) + A2 (X-C1) (X-C2) + ... + an (X-C1) (X-C2) ... (X-CN)

here A0,...,an is the coefficient , while C0,...,CN is the median . Notice that when the C0,...,CN is all 0 o'clock, the Newton interpolation polynomial is degraded to the n-order polynomial defined earlier.

Building an interpolation polynomial

Let's look at how to build an interpolation polynomial on function f (x).

in order to interpolate the function f (x), a polynomial pn (z) with a order less than or equal to n is constructed, which in turn requires a n+1 known point of the function f (x): X0,...,xn. These known point X0,...,xn are called interpolation points . The approximate value of function f (x) at X=z can be calculated by interpolating the polynomial pn (z). the interpolation method needs to satisfy point z within [X0,XN] . The interpolation polynomial pn (z) can be constructed using the following formula.

PN (Z) = f[x0] + f[x0,x1] (z-x0) + f[x0,x1,x2] (z-x0) (z-x1) + ... + f[x0,..., xn] (z-x0) (z-x1) ... (z-xn-1)

where the function f (x) is known at the point x0,..., xn, and f[x0],..., f[x0,..., xn] is called the difference quotient .

The difference quotient can be calculated by point x0,..., xn and the value of function f (x) at these points. This is the calculation formula of Newton interpolation polynomial. Notice the same point of the formula as the Newton interpolation polynomial. The calculation formula for the difference quotient is:

F[xi,..., XJ] = f (xi) if i=j

F[xi,..., XJ] = (f[xi+1,... XJ)-F[xi,... xj-1])/(XJ-XI) if I < J

It is not difficult to see from this formula that when I < J, other difference quotient values must be calculated beforehand. For example, to calculate f[x0,x1,x2,x3], you need to calculate the values of f[x1,x2,x3] and f[x0,x1,x2] first. Fortunately, we can use a bad quotient table to help us calculate the difference quotient value in a systematic way. Such as.

The difference table consists of multiple lines. The top row holds the value of the known point x0,..., xn. The second row holds the value of f[x0],..., F[xn]. to calculate the other difference quotient values in the table, draw a diagonal line from each of the difference quotient values to return to F[XI] and F[XJ] (such as the dashed line at the difference quotient f[x1,x2,x3). to obtain the XI and XJ in the denominator, they are obtained directly from XI and XJ. The two difference quotient in the molecule is the result of the previous stage calculation .

When the calculation of the entire differential quotient table is completed, the coefficients of the interpolation polynomial are the one that starts at the second line and the leftmost item in each row .

Calculating interpolation polynomial

once the coefficients of the interpolation polynomial are determined, for the function f, if we want to know the function value at a point, we only need to evaluate the polynomial .

For example, the function value known as function f at 4 points: x0=-3.0,f (x0) =-5.0;x1=-2.0,f (x1) =-1.1;x2=2.0,f (x2) =1.9;x3=3.0,f (x3) = 4.8; now ask for a point z0=-2.5,z1= The value of the function at the 0.0,z2=1.0,z3=2.5. Since the value of function f at 4 points is known, the interpolation polynomial is 3-step. is a 3-step interpolation polynomial P3 (z) of the difference quotient table.

once the coefficients are obtained from the difference quotient table, the Newton formula described above can be used to construct the interpolation polynomial P3 (z):

P3 (z) =-5.0 + 3.9 (z+3.0) + ( -0.63) (z+3.0) (z+2.0) + 0.1767 (z+3.0) (z+2.0) (z-2.0)

the next step is to calculate the value of the function at each point z by this polynomial. For example, at Point z=-2.5, the following calculations are obtained :

P3 (z) =-5.0 + 3.9 ( -2.5+3.0) + (-0.63) ( -2.5+3.0) ( -2.5+2.0) + 0.1767 ( -2.5+3.0) ( -2.5+2.0) ( -2.5-2.0) = 2.694

The function values at points Z1, Z2, and Z3 can be computed in a similar way. The final result is expressed in tabular and functional images. Such as.

As with any other approximation algorithm, it is important to understand that there are usually some errors associated with the interpolation polynomial. qualitatively, if the error is to be minimized, the built interpolation polynomial must have enough known points on the function f (x) to get to the line. and the distance before the point is appropriate, so that the resulting polynomial can accurately represent the characteristics of the function.

Interface definition for polynomial interpolation

Interpol

int Interpol (const double *x, const double *FX, int n, double *z, double *pz, int m);

return value : If the interpolation operation succeeds, returns 0; otherwise returns-1;

description : The polynomial interpolation method is used to calculate the value of a function at some specific point.

The set of points known by the caller at parameter x specifying the value of the function . the function value corresponding to each known point is specified in FX . The point to be asked is specified by the parameter Z , and the function value corresponding to Z is returned in the PZ . the number of elements in X and F (x) is represented by the parameter n. The number of points to be asked in Z (and the number of returned values in PZ) is represented by the parameter m . The storage space associated with X, FX, Z, and PZ is managed by the caller.

Complexity of: O (Mn2), where m represents the number of values to be evaluated, and N represents the number of known points.

realization and analysis of polynomial interpolation

The polynomial interpolation method is mainly based on the determination of the interpolation polynomial of a series of desired points. To get this polynomial, the coefficients of the polynomial must first be determined by calculating the difference quotient.

first, allocate storage space for the difference quotient and the coefficients to be determined . Note that since each item in each row in the quotient table relies only on the results of its previous row, it is not necessary to keep all the table entries at once. So, just allocate space for the most space-consuming rows , and the row will have n entries .

Next, use the value in FX to initialize the first row of the difference table . This is to prepare for calculating the third row in the difference table. (The first two lines do not need to be evaluated because the entries in both lines are already saved in X and FX).

The last step in initialization is to save the value of fx[0] in coeff[0] because this is the first factor in the interpolation polynomial .

The process of calculating a bad quotient involves a nested loop , in which we calculate the difference quotient in the loop based on the formula described earlier. in the outer loop, K is used to count which row is being computed (excluding the rows represented by X and FX). In the inner loop ,I represents which entry is being computed in the current row . once a line of entries has been computed,the value in Table[0] becomes the next factor in the interpolation polynomial . Therefore, save the value to Coeff[k] . once all the coefficients of the interpolation polynomial are obtained, the values of each target point in Z can be calculated, and the values are then saved in the PZ.

We name this function Interpol, which has a time complexity of O (Mn2), where m represents the number of elements in Z (also the number of values in PZ), and N represents the number of elements in X (also in FX) . The complexity factor, N2, is obtained by changing each iteration in the J Control loop, which has to be multiplied in the current iteration by one more than the previous round. In other words, when j=1, the term needs to do 1 multiplication, when j=2, the term needs to do 2 times multiplication, the process continues until j=n-1, the term needs to do n-1 times multiplication. In fact, this is the sum of the integers of the 1~n-1, and the computed time is t (n) = (n (n+1)/2)-N, multiplied by a fixed period of time. (This is derived from the well-known formula for calculating arithmetic progression). In the large O notation can be simplified to O (N2). The factor m in O (Mn2) is derived from the process of calculating the polynomial interpolation for each point in Z. In the first nested loop, all the difference quotient is computed, and its complexity is O (N2). Therefore, the final complexity has an additional factor m, which has little effect on the actual complexity.

Example: implementation of polynomial interpolation

/*interpol.c*/#include<stdlib.h>#include<string.h>#include"Nummeths.h"/*Interpol*/intInterpolConst Double*x,Const Double*FX,intNDouble*z,Double*pz,intm) {    Doubleterm, *table, *Coeff; inti,j,k; /*allocating space for the difference quotient and the coefficients to be determined*/    if(table = (Double*)malloc(sizeof(Double) (*n)) = =NULL)return-1; if((Coeff = (Double*)malloc(sizeof(Double) (*n)) = =NULL) {         Free(table); return-1; }        /*Initialize Difference table*/memcpy (TABLE,FX,sizeof(Double)*N); /*key points: Determining the coefficients of the difference quotient table*/coeff[0] = table[0];  for(k=1; k<n; k++)    /*the outer loop k is used to count which line is being computed*/    {         for(i=0; i<n-k; i++)        /*The inner Loop i indicates which entry is being computed in the current row (with the increase in the number of rows, the entry is reduced)*/{J=i+K; /*The quotient of molecules in each item of the current line is the result of its previous phase calculation*/Table[i]= (table[i+1]-table[i])/(X[j]-X[i]); }        /*The first entry of the current row evaluates to the next factor of the polynomial*/Coeff[k]=table[0]; }     Free(table); /*evaluates an interpolation polynomial at a specified point (cyclic construction interpolation polynomial)*/     for(k=0; k<m; k++)  /* Outermost: Traverse the Z-point group * /  {/*the first factor of an interpolation polynomial*/Pz[k]= coeff[0];  for(j=1; j<n; J + +)  /* nesting: Constructing a polynomial (new formula equals the result of the previous step plus a new factor)* /  { term=Coeff[j]; /* Factor composition: Based on polynomial coefficients */ for(i=0; i<j; i++)  /* Nested: The result of the new factor above is multiplied by (Z[k]- x[i])* /   term=term* (Z[k]-X[i]); PZ[K]=PZ[K] +term ; }    }         Free(Coeff); return 0;}

Implementation and analysis of polynomial interpolation algorithm based on numerical calculation algorithm