计算方法：拉格朗日插值及其代码实现

Posted on Edited on In , Views:

这篇博文记录计算方法中拉格朗日插值法的Python和matlab实现

数学原理

拉格朗日插值法表达式为：

其中插值基函数

Python实现

LagranggeInterpolation类：

其中定义了fit_interp()和calc_interp()方法分别用于计算拉格朗日多项式、对一个数据点计算插值函数的值

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
import numpy as np
import sympy 

class LagrangeInterpolation:

    def __init__(self, x, y):
        """
        x: 已知数据的x
        y: 已知数据的y
        """

        self.x = np.asarray(x, dtype=float)
        self.y = np.asarray(y, dtype=float)

        if len(self.x) > 1 and len(self.x) == len(self.y):
            self.n = len(self.x)    # 插值点个数
        else:
            return ValueError("插值数据点数量不匹配")
        
        self.polynomial = None  # 最终的插值多项式
        self.poly_cofficient = None # 最终的系数向量，由高次到低次
        self.cofficient_order = None    # 插值多项式次数

    def fit_interp(self):
        """
        插值过程，生成插值多项式
        """
        t = sympy.Symbol("t")
        self.polynomial = 0.0

        for k in range(self.n):
            fun = self.y[k]

            for j in range(0,k):
                fun *= (t-self.x[j])/(self.x[k]-self.x[j])

            for j in range(k+1,self.n):
                fun *= (t-self.x[j])/(self.x[k]-self.x[j])
            
            self.polynomial += fun
        
        self.polynomial = sympy.expand(self.polynomial) # 多项式展开
        polynomial = sympy.Poly(self.polynomial, t) 
        self.poly_cofficient = polynomial.coeffs()
        self.cofficient_order = polynomial.monoms()

    def calc_interp(self, x0):
        """
        计算插值之后的数值
        """
        if self.polynomial == None:
            return ValueError("未生成插值多项式")
        
        x0 = np.asarray(x0, dtype=float)
        n0 = len(x0)
        y_0 = np.zeros(n0)
        t = self.polynomial.free_symbols.pop()
        for i in range(n0):
            y_0[i] = self.polynomial.evalf(subs={t: x0[i]})
        return y_0

举例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
if __name__=='__main__':
    a = np.linspace(0, 2 * np.pi, 10, endpoint=True)
    b = np.sin(a)
    x0 = np.array([np.pi/2, 2.158, 3.58, 4.784])
    lag_interp = LagrangeInterpolation(x=a, y=b)
    lag_interp.fit_interp()
    print("拉格朗日多项式如下：")
    print(lag_interp.polynomial)
    print("拉格朗日插值多项式系数向量和对应阶次")
    print(lag_interp.poly_cofficient)
    print(lag_interp.cofficient_order)
    y0 = lag_interp.calc_interp(x0)
    print("所求插值点的值：", y0)
    print("精确解是     ：", np.sin(x0))

输出

1
2
3
4
5
6
7
拉格朗日多项式如下：
-2.29369601281319e-6*t**9 + 6.48527268907878e-5*t**8 - 0.00061946678451133*t**7 + 0.00167479859899278*t**6 + 0.0040407632285091*t**5 + 0.00707566335693954*t**4 - 0.173843738744177*t**3 + 0.00401585715826069*t**2 + 0.999072579744674*t
拉格朗日插值多项式系数向量和对应阶次
[-2.29369601281319e-6, 6.48527268907878e-5, -0.000619466784511330, 0.00167479859899278, 0.00404076322850910, 0.00707566335693954, -0.173843738744177, 0.00401585715826069, 0.999072579744674]
[(9,), (8,), (7,), (6,), (5,), (4,), (3,), (2,), (1,)]
所求插值点的值： [ 0.99999836  0.83249341 -0.42449807 -0.99743583]
精确解是     ： [ 1.          0.8324932  -0.42449797 -0.99743703]