【问题标题】:How to complete Lagrange Interpolation correctly by pythonpython如何正确完成拉格朗日插值
【发布时间】:2022-12-09 18:32:27
【问题描述】:

我被要求使用拉格朗日插值法绘制一条穿过图形上几个点的线(禁止使用 scipy.interpolate.lagrange)。

但似乎我的代码不能给出正确的结果。说实话,我什至不知道哪一部分是错的。

如果我被要求找到线上特定点的值,我应该如何调整代码来做到这一点?

第一张图片是预期的输出。第二个是我的代码的输出。

有人可以告诉我如何纠正吗?

源代码

import numpy as np
import matplotlib.pyplot as plt

m = 4
X = np.array([0, 1, 1.5, 2, 3])
Y = np.array([2, 8, 12, 10, 20])

def p(x): 
    px = 0;
    for j in range(m+1):
        Lmj = 1
        for k in range(m+1):
            if k != j:
                Lmj *= (x - X[k])/(X[j] - X[k])
            px += Y[j] * Lmj;
        return px

plt.scatter(X, Y)
Xinterp = np.linspace(min(X), max(X), 100)
plt.plot(Xinterp, p(Xinterp))
plt.show()

另外, plt.show() 行是否必要?从我讲师的笔记来看,代码没有这一行,但仍然可以显示图形。但是在我复制的源代码版本中,我需要添加这一行来显示图形。

非常感谢您注意到这个问题。

【问题讨论】:

    标签: python interpolation graphplot


    【解决方案1】:
    1. 拉格朗日插值从一组给定的点生成多项式。
    2. 然后您可以使用此多项式求出直线上特定点的值。
    3. 如果您在 Jupyter notebook(或类似笔记本)中运行代码,则 plt.show() 不是必需的,但如果通过终端运行则需要 plt.show() 。

      以下是一些 python 代码,可根据您的点集生成 numpy.polynomial

      import numpy as np
      from numpy.polynomial import Polynomial
      import matplotlib.pyplot as plt
      
      X = np.array([0, 1, 1.5, 2, 3])
      Y = np.array([2, 8, 12, 10, 20])
      
      n = len(X)
      poly = Polynomial(np.zeros(n))
      
      for j in range(n):
          k = [k for k in range(n) if k != j]
          roots = -1 * X[k]
      
          sub_poly = Polynomial.fromroots(X[k])
          scale = Y[j] / np.prod(X[j] - X[k])
          sub_poly.coef *= scale
      
          poly.coef += sub_poly.coef
      

      然后您可以通过以下方式绘制它:

      plt.scatter(X, Y)
      Xinterp = np.linspace(min(X), max(X), 100)
      plt.plot(Xinterp, poly(Xinterp))
      plt.show()
      

      或查找特定值:

      print(poly(0.5))
      

      更多信息请看这里:https://learn.64bitdragon.com/articles/computer-science/numerical-analysis/lagrange-interpolation

      完全披露:我是这篇文章的作者。

    【讨论】:

      猜你喜欢
      • 2015-01-20
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2019-09-23
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      相关资源
      最近更新 更多