2021-09-15

python数值分析(一、拉格朗日插值法)

最近在看论文中，发现地信和遥感的很多算法都脱离不开数值分析那一套，所以我认为要想更好的处理以后数据和工作，在数学模型算法实现上必须狠下功夫，希望我能利用寒假这段时间，把数值分析中一些算法用python实现，以此谨记!

# -*- coding: utf-8 -*-
"""
# @Time    : 2022/1/15 18:59
# @Author  : FJC
# @File    : test_lagrange.py
# @Software: win10  python3.7
"""

import numpy as np
import matplotlib.pyplot as plt
import sympy


class LagrangeInterpolation:
    """
    拉格朗日插值法
    """

    def __init__(self, x, y):
        """
        初始化变量
        :param x:已知数据的x坐标
        :param y:已知数据的y坐标
        """
        self.x = x
        self.y = y
        self.n = len(self.x)  # 有多少已知点就是几次基函数
        self.polynomial = None

    def lagrange_interpolation(self):
        """
        拉格朗日插值法核心算法
        :return:
        """
        t = sympy.Symbol("t")
        self.polynomial = 0.0
        # 三个点是两次基函数，即两次基函数要三个点
        for i in range(self.n):
            basis_fun = self.y[i]  # 将已知y坐标赋值
            for j in range(i):
                basis_fun *= (t - self.x[j]) / (self.x[i] - self.x[j])
            for j in range(i + 1, self.n):
                basis_fun *= (t - self.x[j]) / (self.x[i] - self.x[j])
            self.polynomial += basis_fun
        print(self.polynomial)

    def cal_interp_x0(self, x1):
        """
        插值未知x坐标
        :param x1:需要插值的数据
        :return:
        """
        x1 = np.asarray(x1, dtype=np.float64)
        n = len(x1)
        y1 = np.zeros(n)
        t = self.polynomial.free_symbols
        t = self.polynomial.free_symbols.pop()
        for i in range(n):
            y1[i] = self.polynomial.subs({t: x1[i]})
        return y1

    def plt_lagrange(self, x2):
        """
        画图
        :param x2: 插值100个数据
        :return:
        """
        y2 = self.cal_interp_x0(x2)
        plt.figure(figsize=(8, 6))
        plt.plot(self.x, self.y, "ro", label="Interpolation base points")
        plt.plot(x2, y2, "b--", label="Interpolation polynomial")
        plt.legend()
        plt.xlabel("x", fontdict={"fontsize": 12})
        plt.ylabel("y", fontdict={"fontsize": 12})
        plt.title("Lagerange interpolation", fontdict={"fontsize": 14})
        plt.show()


if __name__ == "__main__":
    x = np.linspace(0, 2 * np.pi, 10)
    y = np.sin(x)
    x1 = np.array([2.456, 3.254, 6.345])  # 测试需要插值的数据
    x2 = np.linspace(0, 2 * np.pi, 100)
    print(y)
    lagrange = LagrangeInterpolation(x, y)
    lagrange.lagrange_interpolation()
    print(lagrange.cal_interp_x0(x1))
    lagrange.plt_lagrange(x2)