Python实现数值计算----牛顿插值法

  拉格朗日插值法的最大毛病就是每次引入一个新的插值节点,基函数都要发生变化,这在一些实际生产环境中是不合适的,有时候会不断的有新的测量数据加入插值节点集,

因此,通过寻找n个插值节点构造的的插值函数与n+1个插值节点构造的插值函数之间的关系,形成了牛顿插值法。推演牛顿插值法的方式是归纳法,也就是计算Ln(x)- Ln+1(x),并且从n=1开始不断的迭代来计算n+1时的插值函数。

  牛顿插值法的公式是:

  注意:在程序中我用W 代替 

  计算牛顿插值函数关键是要计算差商,n阶差商的表示方式如下:

                        

    关于差商我在这里并不讨论

  计算n阶差商的公式是这样:

  很明显计算n阶差商需要利用到两个n-1阶差商,这样在编程的时候很容易想到利用递归来实现计算n阶差商,不过需要注意的是递归有栈溢出的潜在危险,在计算差商的时候

更是如此,每一层递归都会包含两个递归,递归的总次数呈满二叉树分布:

    

  这意味着递归次数会急剧增加:(。所以在具体的应用中还需要根据应用来改变思路或者优化代码

  废话少说,放码过来。

  首先写最关键的一步,也就是计算n阶差商:

"""
@brief:   计算n阶差商 f[x0, x1, x2 ... xn]
@param:   xi   所有插值节点的横坐标集合                                                        o
@param:   fi   所有插值节点的纵坐标集合                                                      /   @return:  返回xi的i阶差商(i为xi长度减1)                                                     o     o
@notice:  a. 必须确保xi与fi长度相等                                                        / \   /           b. 由于用到了递归,所以留意不要爆栈了.                                           o   o o   o
          c. 递归减递归(每层递归包含两个递归函数), 每层递归次数呈二次幂增长,总次数是一个满二叉树的所有节点数量(所以极易栈溢出)
"""
def get_order_diff_quot(xi = [], fi = []):
    if len(xi) > 2 and len(fi) > 2:
        return (get_order_diff_quot(xi[:len(xi) - 1], fi[:len(fi) - 1]) - get_order_diff_quot(xi[1:len(xi)], fi[1:len(fi)])) / float(xi[0] - xi[-1])
    return (fi[0] - fi[1]) / float(xi[0] - xi[1])

  看上面的牛顿插值函数公式,有了差商,还差

  这个就比较好实现了:

"""
@brief:  获得Wi(x)函数;
         Wi的含义举例 W1 = (x - x0); W2 = (x - x0)(x - x1); W3 = (x - x0)(x - x1)(x - x2)
@param:  i  i阶(i次多项式)
@param:  xi  所有插值节点的横坐标集合
@return: 返回Wi(x)函数
"""
def get_Wi(i = 0, xi = []):
    def Wi(x):
        result = 1.0
        for each in range(i):
            result *= (x - xi[each])
        return result
    return Wi

    

    OK, 牛顿插值法最重要的两部分都有了,下面就是将这两部分组合成牛顿插值函数,如果是c之类的语言就需要保存一些中间数据了,我利用了Python的闭包直接返回一个牛顿插值函数,闭包可以利用到它所处的函数之中的上下文数据。

"""
@brief: 获得牛顿插值函数
@
"""
def get_Newton_inter(xi = [], fi = []):
    def Newton_inter(x):
        result = fi[0]
        for i in range(2, len(xi)):
            result += (get_order_diff_quot(xi[:i], fi[:i]) * get_Wi(i-1, xi)(x))
        return result
    return Newton_inter

    上面这段代码就是对牛顿插值函数公式的翻译,注意get_Wi函数的参数是i-1,这个从函数的表达式可以找到原因。

  构造一些插值节点

 

 ‘‘‘ 插值节点, 这里用二次函数生成插值节点,每两个节点x轴距离位10 ‘‘‘
    sr_x = [i for i in range(-50, 51, 10)]
    sr_fx = [i**2 for i in sr_x]

   

 获得牛顿插值函数

    Nx = get_Newton_inter(sr_x, sr_fx)            # 获得插值函数

    tmp_x = [i for i in range(-50, 51)]          # 测试用例
    tmp_y = [Nx(i) for i in tmp_x]               # 根据插值函数获得测试用例的纵坐标

  

完整代码:

# -*- coding: utf-8 -*-
"""
Created on Thu Nov 17 18:34:21 2016

@author: tete
@brief: 牛顿插值法
"""

import matplotlib.pyplot as plt

"""
@brief:   计算n阶差商 f[x0, x1, x2 ... xn]
@param:   xi   所有插值节点的横坐标集合                                                        o
@param:   fi   所有插值节点的纵坐标集合                                                      /   @return:  返回xi的i阶差商(i为xi长度减1)                                                     o     o
@notice:  a. 必须确保xi与fi长度相等                                                        / \   /           b. 由于用到了递归,所以留意不要爆栈了.                                           o   o o   o
          c. 递归减递归(每层递归包含两个递归函数), 每层递归次数呈二次幂增长,总次数是一个满二叉树的所有节点数量(所以极易栈溢出)
"""
def get_order_diff_quot(xi = [], fi = []):
    if len(xi) > 2 and len(fi) > 2:
        return (get_order_diff_quot(xi[:len(xi) - 1], fi[:len(fi) - 1]) - get_order_diff_quot(xi[1:len(xi)], fi[1:len(fi)])) / float(xi[0] - xi[-1])
    return (fi[0] - fi[1]) / float(xi[0] - xi[1])

"""
@brief:  获得Wi(x)函数;
         Wi的含义举例 W1 = (x - x0); W2 = (x - x0)(x - x1); W3 = (x - x0)(x - x1)(x - x2)
@param:  i  i阶(i次多项式)
@param:  xi  所有插值节点的横坐标集合
@return: 返回Wi(x)函数
"""
def get_Wi(i = 0, xi = []):
    def Wi(x):
        result = 1.0
        for each in range(i):
            result *= (x - xi[each])
        return result
    return Wi

"""
@brief: 获得牛顿插值函数
@
"""
def get_Newton_inter(xi = [], fi = []):
    def Newton_inter(x):
        result = fi[0]
        for i in range(2, len(xi)):
            result += (get_order_diff_quot(xi[:i], fi[:i]) * get_Wi(i-1, xi)(x))
        return result
    return Newton_inter

"""
demo:
"""
if __name__ == ‘__main__‘:   

    ‘‘‘ 插值节点, 这里用二次函数生成插值节点,每两个节点x轴距离位10 ‘‘‘
    sr_x = [i for i in range(-50, 51, 10)]
    sr_fx = [i**2 for i in sr_x]  

    Nx = get_Newton_inter(sr_x, sr_fx)            # 获得插值函数

    tmp_x = [i for i in range(-50, 51)]          # 测试用例
    tmp_y = [Nx(i) for i in tmp_x]               # 根据插值函数获得测试用例的纵坐标

    ‘‘‘ 画图 ‘‘‘
    plt.figure("I love china")
    ax1 = plt.subplot(111)
    plt.sca(ax1)
    plt.plot(sr_x, sr_fx, linestyle = ‘‘, marker=‘o‘, color=‘b‘)
    plt.plot(tmp_x, tmp_y, linestyle = ‘--‘, color=‘r‘)
    plt.show()

  

  

时间: 2024-10-12 08:56:08

Python实现数值计算----牛顿插值法的相关文章

python与数值计算环境搭建

数值计算的编程的软件很多种,也见过一些编程绘图软件的对比. 利用Python进行数值计算,需要用到numpy(矩阵) ,scipy(公式符号), matplotlib(绘图)这些工具包. 1.Linux系统中一般会带有Python.可以用命令查看是否安装Python $ python Python 2.7.5 (default, Feb 11 2014, 07:46:25) [GCC 4.8.2 20140120 (Red Hat 4.8.2-13)] on linux2 Type "help&

牛顿插值法及其C++实现

h1 { margin-bottom: 0.21cm } h1.western { font-family: "Liberation Sans", sans-serif; font-size: 18pt } h1.cjk { font-family: "Noto Sans CJK SC Regular"; font-size: 18pt } h1.ctl { font-family: "Noto Sans CJK SC Regular"; fon

Python实现牛顿插值法(差商表)

def func(x,y,X,infor=True): list2=[y[0]] # 差商表的对角线的第一个元素始终是y0 count=1 while(True): if len(y)>1: list=[] # 空列表用来保存,每次计算后差商表的行 for i in range(len(y)-1): n=x[i+count]-x[i] m=y[i+1]-y[i] l=m/n list.append(l) list2.append(list[0]) # list2用来记录差商表的对角线元素,每计算

复化梯形求积分——用Python进行数值计算

用程序来求积分的方法有很多,这篇文章主要是有关牛顿-科特斯公式. 学过插值算法的同学最容易想到的就是用插值函数代替被积分函数来求积分,但实际上在大部分场景下这是行不通的. 插值函数一般是一个不超过n次的多项式,如果用插值函数来求积分的话,就会引进高次多项式求积分的问题.这样会将原来的求积分问题带到另一个求积分问题:如何求n次多项式的积分,而且当次数变高时,会出现龙悲歌现象,误差反而可能会增大,并且高次的插值求积公式有可能会变得不稳定:详细原因不赘述. 牛顿-科特斯公式解决这一问题的办法是将大的插

Python 3 数值计算

Python 3.4.3 (v3.4.3:9b73f1c3e601, Feb 24 2015, 22:43:06) [MSC v.1600 32 bit (Intel)] on win32Type "copyright", "credits" or "license()" for more information.>>> 17 /3 #典型的除法返回一个浮点数5.666666666666667>>> 17 //

用python做数值计算

http://sebug.net/paper/books/scipydoc/scipy_intro.html http://www.cnblogs.com/weilq/p/3432817.html https://eclipse.org/downloads/packages/eclipse-classic-422/junosr2 https://www.ics.uci.edu/~pattis/common/handouts/introtopythonineclipse/

Python实现数值计算----分段二次插值

事实上在实际使用中,高次插值显然是很不适合的,高次插值将所有样点包涵进一个插值函数中,这是次幂高的原因.高次计算复杂,而且刚开始的一点误差会被方的很大.因此将整个区间分为若干个小区间,在每一个小区间进行插值这样更好,实现容易,也方便在一些嵌入式设备上使用.有不少需要插值方法的场景是在嵌入式的应用中. 我以等距节点的二次插值为例,以每三个节点为一个子区间. 等距节点二次插值很好写,由于每个区间只有三个插值节点,计算差商也不必使用拉格朗日插值中使用的递归,直接列表达式也很简单(实际上等距节点二次插值

【数值分析】拉格朗日插值与牛顿插值

在工程应用和科学研究中,经常要研究变量之间的关系y=f(x).但对于函数f(x),常常得不到一个具体的解析表达式,它可能是通过观测或实验得到的一组数据(x,f(x)),x为一向量;或则是解析表达式非常复杂,不便于计算和使用.因此我们需要寻找一个计算比较简单的函数S(x)近似代替f(x),并使得S(x)=f(x),这种方法就称为插值法. 常用的插值法有: 一维插值法:拉格朗日插值.牛顿插值.分段低次插值.埃尔米特插值.样条插值. 二维插值法:双线性插值.双二次插值. 拉格朗日插值法 已知函数f(x

拉格朗日插值法(图文详解)

在数值分析中,拉格朗日插值法是以法国十八世纪数学家约瑟夫·拉格朗日命名的一种多项式插值方法.许多实际问题中都用函数来表示某种内在联系或规律,而不少函数都只能通过实验和观测来了解.如对实践中的某个物理量进行观测,在若干个不同的地方得到相应的观测值,拉格朗日插值法可以找到一个多项式,其恰好在各个观测的点取到观测到的值.这样的多项式称为拉格朗日(插值)多项式.数学上来说,拉格朗日插值法可以给出一个恰好穿过二维平面上若干个已知点的多项式函数.拉格朗日插值法最早被英国数学家爱德华·华林于1779年发现[1