[Numerical calculation methods] Curve fitting and interpolation: Lagrange interpolation, Newton interpolation and their python/C implementation

Table of Contents

1. Approximate expression

Interpolation

Fitting

Projection

2. Interpolation

1. Lagrange interpolation

Lagrange interpolation formula

Linear interpolation(n=1)

Parabolic interpolation (n=2)

python implementation

C language implementation

2. Newton interpolation

python implementation

C language implementation

1. Approximate expression

Interpolation, fitting, and projection are common approximate expressions used to estimate, predict, or represent data or functions.

  • Interpolation

    • Refers to estimating or deducing the values between known data points through interpolation between these data points. Interpolation can be used to construct smooth curves or surfaces to make predictions between data points or to supplement missing data.

  • Fitting

    • It refers to the process of best fitting a mathematical model to known data points by selecting appropriate functional forms and parameters. The goal of fitting is to find a functionwhose values near the data points are as close as possible to the actual observed values. Fitting can be used in data analysis, curve fitting, regression analysis and other fields.

  • Projection

    • Refers to the process of mapping a vector or a set of vectors to another vector space or subspace. In linear algebra, projection can be used to find the projection or projected components of a vector onto another vector or vector space. Projection can be used in fields such as dimensionality reduction, data compression, and feature extraction, as well as projection transformation in computer graphics.

2. Interpolation

Lagrange interpolation and Newton interpolation are common polynomial interpolation methods used to estimate function values at other points through a given set of data points. The main difference between them is how the interpolation polynomial is constructed.

  • Lagrange interpolation uses a method based on Lagrange polynomials to construct interpolation polynomials.
    • Lagrange polynomials are constructed by multiplying each data point by a basis function such that the basis function is zero at other data points. The final interpolation polynomial is obtained by adding all these basis functions.
    • The advantage of Lagrange interpolation is that it is easy to understand and implement, but it may cause high computational complexity when there are many data points.
  • Newton interpolation uses the concept of difference quotient to construct the interpolation polynomial.
    • The difference quotient is a recursively defined concept used to calculate the coefficients in interpolating polynomials. The calculation of the difference quotient can be performed in tabular form, where each column represents a different order of difference quotients. By calculating the difference quotient, the interpolation polynomial can be constructed step by step.
    • The advantage of Newton interpolation is that the calculated difference quotient value can be reused when calculating the difference quotient, thereby reducing the amount of calculation.

1. Lagrange interpolation

Lagrange interpolation is a method for constructing a polynomial function from known data points. It is based on the principle of Lagrangian interpolation polynomial, which passes through each data point and satisfies the corresponding conditions. Lagrangian interpolation can be used to estimate values between data points rather than just interpolating at a given data point.

The basic steps for using Lagrange interpolation are as follows:

  • Given a set of known data points, including the values of the abscissa and ordinate.
  • According to the number of data points, a Lagrangian interpolation polynomial of corresponding degree is constructed.
  • Multiply the function value of each data point by the corresponding Lagrangian interpolation polynomial and add them to get the final interpolation function.

In this way, a smooth interpolation function can be obtained at a given data point, such that the value of the function can be estimated at any position between these data points. Lagrange interpolation may have some problems when there are few data points or there are large gaps between data points. For example, the interpolation polynomial may produce oscillation, which is called the Runge phenomenon.

Lagrange interpolation formula

Linear interpolation (n=1)

Parabolic interpolation (n=2)

Vandermonde determinant – Zhihu (zhihu.com)icon-default.png?t=N7T8https://zhuanlan.zhihu.com/p/161300510

python implementation

import numpy as np


#Define Lagrange interpolation function
def lagrange_interpolation(x, y, xi):
    n = len(x)
    yi = 0.0

    for i in range(n):
        # Calculate each term of the Lagrangian interpolation polynomial
        term = y[i]
        for j in range(n):
            if j != i:
                term *= (xi - x[j]) / (x[i] - x[j])
        yi + = term

    return yi


# Example data points
x = np.array([0.32, 0.34, 0.36])
y = np.array([0.314567, 0.333487, 0.352274])

# Points to be interpolated
xi = 0.3367

# Perform interpolation
yi = lagrange_interpolation(x, y, xi)

print("Interpolation result:", yi)
print("Real result:", np.sin(xi))

Output:

Interpolation result: 0.3303743620374999
Real result: 0.330374191555628

C language implementation

#include <stdio.h>

// Calculate the value of Lagrange interpolation polynomial
double lagrange_interpolation(double x[], double y[], int n, double xi) {
    double yi = 0.0;

    for (int i = 0; i < n; i + + ) {
        double term = y[i];
        for (int j = 0; j < n; j + + ) {
            if (j != i) {
                term *= (xi - x[j]) / (x[i] - x[j]);
            }
        }
        yi + = term;
    }

    return yi;
}

int main() {
    //Example data points

    double x[] = {0.32, 0.34, 0.36};
    double y[] = {0.314567, 0.333487, 0.352274};

    //Points to be interpolated
    double xi = 0.3367;

    //The number of data points
    int n = sizeof(x) / sizeof(x[0]);

    // perform interpolation
    double yi = lagrange_interpolation(x, y, n, xi);

    printf("Interpolation result: %f\\
", yi);

    return 0;
}

Output:

Interpolation result: 0.330374

2. Newton interpolation

Newton interpolation is based on the concept of difference quotient. Given a set of data points, Newton interpolation can generate a polynomial passing through these points to perform interpolation and extrapolation within the given data range.

The basic idea of Newton interpolation is to use the difference quotient to recursively construct a polynomial. The difference quotient is defined by recursively calculating the difference between data points. Specifically, for the given data points (x0, y0), (x1, y1), …, (xn, yn), the difference quotient can be expressed as:

f[x_{0}] = y_{0}
f[x_{1}, x_{0}] =\frac{ (f[x_{1}] - f[x_{0}]) }{ (x_{1} - x_{ 0})}

f[x_{2}, x_{1}, x_{0}] =\frac{ (f[x_{2}, x_{1}] - f[x_{1}, x_ {0}]) }{ (x_{2} - x_{0})}
…………
f[xn, xn-1, …, x0] = (f[xn, xn-1, …, x1] – f[xn-1, …, x0]) / (xn – x0)

Then, by adding these difference quotients step-by-step to the polynomial, you get a polynomial expressed as:

P(x) = f[x_{0}] + f[x_{1}, x_{0}](x - x_{0}) + f[x_{2}, x_{1 }, x_{0}](x - x_{0})(x - _{x1}) + ...

This polynomial can be used to interpolate within a given range of data, i.e. using known data points to estimate the value of a function at other points. One of the advantages of Newton interpolation is that it can incrementally improve the interpolation results by adding more data points. However, it also has some problems. For example, the resulting polynomial may suffer from Runge’s phenomenon, causing oscillations at the boundaries.

python implementation

def newton_interpolation(x, y, xi):
    # Calculate the difference quotient
    n = len(x)
    f = [[0] * n for _ in range(n)]
    for i in range(n):
        f[i][0] = y[i]

    for j in range(1, n):
        for i in range(n - j):
            f[i][j] = (f[i + 1][j - 1] - f[i][j - 1]) / (x[i + j] - x[i])

    # Construct interpolation polynomial
    result = f[0][0]
    for j in range(1, n):
        term = f[0][j]
        for i in range(j):
            term *= (xi - x[i])
        result + = term

    return result


# Sample data
x = [0.32, 0.34, 0.36]
y = [0.314567, 0.333487, 0.352274]
xi = 0.3367

# Perform interpolation
interpolated_value = newton_interpolation(x, y, xi)
print("Interpolation result:", interpolated_value)

Output:

Interpolation result: 0.3303743620375

C language implementation

#include <stdio.h>

double newton_interpolation(double x[], double y[], int n, double xi) {
    // Calculate the difference quotient
    double f[n][n];
    for (int i = 0; i < n; i + + ) {
        f[i][0] = y[i];
    }

    for (int j = 1; j < n; j + + ) {
        for (int i = 0; i < n - j; i + + ) {
            f[i][j] = (f[i + 1][j-1] - f[i][j-1]) / (x[i + j] - x[i]);
        }
    }

    // Construct interpolation polynomial
    double result = f[0][0];
    for (int j = 1; j < n; j + + ) {
        double term = f[0][j];
        for (int i = 0; i < j; i + + ) {
            term *= (xi - x[i]);
        }
        result + = term;
    }

    return result;
}

int main() {
    // sample data
    double x[] = {0.32, 0.34, 0.36};
    double y[] = {0.314567, 0.333487, 0.352274};
    int n = sizeof(x) / sizeof(x[0]);
    double xi = 0.3367;

    // perform interpolation
    double interpolated_value = newton_interpolation(x, y, n, xi);
    printf("Interpolation result: %f\\
", interpolated_value);

    return 0;
}

Output:

Interpolation result: 0.330374