在文章拉格朗日插值多项式的原理介绍及其应用中,笔者介绍了如何使用拉格朗日插值多项式来拟合任意数据点集。
事实上,插值多项式会更倾向于某些形状。德国数学家卡尔·龙格Carl Runge发现,插值多项式在差值区间的端点附近会发生扭动,且波动较大。这就是数值分析中著名的龙格现象(Runge Phenomenon)。
本文以函数
f
(
x
)
=
1
1
+
12
x
2
f(x)=\frac{1}{1+12x^{2}}
f(x)=1+12x21和区间[-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越大,扭动的幅度就越大,这是用计算机程序对龙格现象的一个模拟。