|
- # !/usr/bin/env python
- # -*- coding: utf-8 -*-
- """
- @Time : 2020/12/30 19:32
- @Author : Albert Darren
- @Contact : 2563491540@qq.com
- @File : my_interpolation.py
- @Version : Version 1.0.0
- @Description : TODO 自己实现拉格朗日插值,牛顿均差插值,埃尔米特插值,分段低次插值,三次样条插值
- @Created By : PyCharm
- """
- import numpy as np
- from sympy import expand, Poly, cos
- from sympy.abc import x
-
-
- # expr = x*y + x - 3 + 2*x**2 - z*x**2 + x**3
- # print(collect(expr, x)) #按x降幂排列
- # print(collect(expr, x).coeff(x,2)) #和 coeff搭配使用,可以求n次项系数, 这里就是求x^2项的系数
-
-
- # 自己原创
- # 拉格朗日插值多项式输出显示形式与scipy.interpolate模块的同名函数有差距,有待继续优化
- def lagrange(x_list: list, y_array: np.ndarray):
- """
- 实现拉格朗日插值函数
- :param x_list: 插值点横坐标列表
- :param y_array: 插值点纵坐标数组
- :return: 拉格朗日插值多项式
- """
- degree = y_array.shape[0]
- l_n = 0
- for k in range(degree):
- # 注意不能同时对列表进行遍历和删除,否则索引会超出范围报错,有几种方式解决,此处使用切片
- l_n += y_array[k] * lagrange_base_function(k, x_list[::])
- return expand(l_n)
-
-
- # 自己原创
- def lagrange_base_function(k, x_list: list):
- """
- 实现拉格朗日插值基函数
- :param k: 插值结点横坐标x_k的下标
- :param x_list: 插值点横坐标列表
- :return: 第k个插值结点的插值基函数
- """
- x_k = x_list[k] # 取得xk处的值
- x_list.pop(k) # 删除xk处的值
- x_list = np.array(x_list)
- # 向量运算实现插值基函数
- return np.prod((x - x_list) / (x_k - x_list))
- # 按元素实现插值基函数
- # l_k = 1
- # for i in range(x_coefficient.shape[0]):
- # l_k *= (x - x_coefficient[i]) / (x_k - x_coefficient[i])
- # return l_k
-
-
- # 自己原创
- def newtown_base_function(symbol_x, x_array: np.ndarray):
- """
- 实现牛顿均差插值基函数
- :param symbol_x: 符号变量x
- :param x_array: 横坐标数组
- :return: 牛顿均差插值基函数数组
- """
- nth = x_array.shape[0]
- base_function_array = np.zeros(nth, dtype=np.ndarray)
- # base_function_array = list(np.zeros(nth))
- base_function_array[0] = 1
- for i in range(1, nth):
- base_function_array[i] = np.prod(symbol_x - x_array[:i], dtype=np.ndarray)
- return base_function_array
-
-
- # 自己原创
- def difference_quotient_coefficient(x_array: np.ndarray, y_array: np.ndarray):
- """
- 实现牛顿均差插值多项式各阶均差系数
- :param x_array: 横坐标数组
- :param y_array: 纵坐标数组
- :return: 各阶均差系数数组
- """
- nth = x_array.shape[0]
- # 初始化均差系数数组
- difference_quotient_array = np.zeros(nth)
- difference_quotient_array[0] = y_array[0]
- # 求各阶均差
- for i in range(1, nth):
- nth_difference_quotient = np.zeros(nth - i)
- # 求每阶均差
- for j in range(nth - i):
- nth_difference_quotient[j] = (y_array[j + 1] - y_array[j]) / (x_array[i + j] - x_array[j])
- difference_quotient_array[i] = nth_difference_quotient[0]
- y_array = nth_difference_quotient
- return difference_quotient_array
-
-
- # 自己原创
- def newtown(x_array: np.ndarray, y_array: np.ndarray, symbol_x):
- """
- 实现牛顿均差插值多项式
- :param x_array: 插值结点横坐标数组
- :param y_array: 插值结点纵坐标数组
- :param symbol_x: 符号变量x
- :return: 牛顿均差插值多项式
- """
- return expand(np.sum(difference_quotient_coefficient(x_array, y_array) * newtown_base_function(symbol_x, x_array)))
-
-
- # 自己原创
- def qin_jiu_shao(coefficient_array: np.ndarray, value, is_ascending_order: bool = False):
- """
- 多项式求值的秦九韶算法,即Horner算法
- :param coefficient_array: 多项式系数数组
- :param value: 自变量的值
- :param is_ascending_order: 多项式排列顺序
- :return: the value of polynomial
- example:
- 2*x**4 - 3*x**2 + 3*x - 4=((((2x+0)x-3)x+3)x-4)
- 3*x**5 - 2*x**3 + x + 7=(((((3x+0)x-2)x+0)x+1)x+7)
- """
- if is_ascending_order:
- coefficient_array = coefficient_array[::-1]
- degree = coefficient_array.shape[0]
- b = coefficient_array[0]
- for i in range(1, degree):
- b = b * value + coefficient_array[i]
- return b
-
-
- if __name__ == '__main__':
- # 拉格朗日插值多项式测试成功,来源详见李庆扬数值分析第5版P48,e.x.1
- """
- x_coefficient_list = [1, -1, 2]
- x_coefficient_array = np.array(x_coefficient_list)
- y = np.array([0, -3, 4])
- lagrange_interpolation_polynomial = lagrange(x_coefficient_list, y)
- # 输出拉格朗日插值多项式
- print("拉格朗日插值多项式为:{}".format(lagrange_interpolation_polynomial))
- # 计算拉格朗日插值多项式在x=2处的值
- print("拉格朗日插值多项式在x=2处的值为:{}".format(lagrange_interpolation_polynomial.subs(x, 2)))
- """
- # 牛顿均差插值多项式测试成功,来源详见李庆扬数值分析第5版P48,e.x.1
- """
- x_coefficient_list = [1, -1, 2]
- x_coefficient_array = np.array(x_coefficient_list)
- y = np.array([0, -3, 4])
- # 牛顿均差插值多项式均差系数测试成功
- print("牛顿均差插值多项式均差系数:{}".format(difference_quotient_coefficient(x_coefficient_array, y)))
- # 牛顿均差插值多项式基函数测试成功
- print("牛顿均差插值多项式基函数:{}".format(newtown_base_function(x, x_coefficient_array)))
- newtown_polynomial = newtown(x_coefficient_array, y, x)
- # 计算牛顿插值多项式(降幂排列)在x=2处的值
- print("牛顿插值多项式(降幂排列)在x=2处的值为:{}".format(newtown_polynomial.subs(x, 2)))
- # 创建牛顿插值多项式对象
- newtown_interpolation_polynomial = Poly(newtown_polynomial, x)
- # 获得牛顿插值多项式降幂排列系数
- newtown_polynomial_coefficient = np.array(newtown_interpolation_polynomial.coeffs()) # list 类型
- # 使用秦九韶算法计算牛顿插值多项式(降幂排列)在x=2处的值
- print("秦九韶算法计算牛顿插值多项式(降幂排列)在x=2处的值为:{}".format(qin_jiu_shao(newtown_polynomial_coefficient, 2)))
- # 牛顿插值多项式次数
- print("牛顿插值多项式次数为:{}".format(newtown_interpolation_polynomial.degree()))
- # 牛顿插值多项式子项次数
- print("牛顿插值多项式子项次数为:{}".format(newtown_interpolation_polynomial.monoms()))
- """
|