插值多项式的龙格现象的介绍与模拟

  在文章拉格朗日插值多项式的原理介绍及其应用中，笔者介绍了如何使用拉格朗日插值多项式来拟合任意数据点集。
  事实上，插值多项式会更倾向于某些形状。德国数学家卡尔·龙格Carl Runge发现，插值多项式在差值区间的端点附近会发生扭动，且波动较大。这就是数值分析中著名的龙格现象（Runge Phenomenon）。
  本文以函数$f(x)=\frac{1}{1+12x^2}$和区间[-1,1]为例，在该区间上平均取n个点（包括端点），在函数图像上得到n个样本点，对这些样本点使用拉格朗日插值多项式，并绘制该插值多项式的图像，观察其在端点附近的表现。
  Python实现程序如下：

# -*- coding: utf-8 -*-
# @Time : 2023/3/8 18:55
# @Author : Jclian91
# @File : runge_phenomenon.py
# @Place : Xuhui, Shanghai
import matplotlib.pyplot as plt


# sample function
# 函数f(x)=1/(1+12*x**2)
def sample_func(x):
    return 1 / (1 + 12 * x ** 2)


# get sample points from sample function with interval [-1, 1]
def get_sample_points(n):
    # n: number of sample points
    step = 2 / (n-1)
    x_values = [-1 + i * step for i in range(n)]
    y_values = [sample_func(x) for x in x_values]
    return x_values, y_values


# get basic lagrange polynomial unit
def get_lagrange_polynomial_unit(x_values, k, x):
    # x_values: values of x in list x_values
    # k: kth lagrange polynomial unit
    # x: variable in kth lagrange polynomial unit
    poly_unit = 1
    for i in range(len(x_values)):
        if i != k:
            poly_unit *= (x-x_values[i])/(x_values[k]-x_values[i])
    return poly_unit


# get lagrange polynomial
def get_lagrange_polynomial(x_values, y_values, x):
    poly = 0
    for i, y in enumerate(y_values):
        poly += y * get_lagrange_polynomial_unit(x_values, i, x)
    return poly


# plot curves with matplotlib
def plot_function(n):
    # plot lagrange polynomial with n sample points from sample function
    sample_x_values, sample_y_values = get_sample_points(n)
    sample_points_number = 500
    x_list = [-1 + i * 2 / (sample_points_number-1) for i in range(sample_points_number)]
    original_y_list = [sample_func(x) for x in x_list]
    y_list = [get_lagrange_polynomial(sample_x_values, sample_y_values, x)for x in x_list]
    plt.plot(x_list, original_y_list, label='f(x)=1/(1+12*x**2)')
    plt.plot(x_list, y_list, label='lagrange polynomial')
    plt.title(f'Runge phenomenon with {n} basic points in function f(x)=1/(1+12*x**2)')
    plt.legend()
    # plt.show()
    plt.savefig(f"{n}_basic_points.png")


if __name__ == '__main__':
    n_points = 5
    plot_function(n_points)

当n=5时，拉格朗日插值多项式的图像如下：

当n=15，拉格朗日插值多项式的图像如下：

当n=25时，拉格朗日插值多项式的图像如下：

当n=35时，拉格朗日插值多项式的图像如下：

当n=45时，拉格朗日插值多项式的图像如下：

通过上述程序的模拟结果，我们可以发现该插值多项式在区间端点附近会发生扭动，当n越大，扭动的幅度就越大，这是用计算机程序对龙格现象的一个模拟。