牛顿插值法——用Python进行数值计算
拉格朗日插值法的最大毛病就是每次引入一个新的插值节点,基函数都要发生变化,这在一些实际生产环境中是不合适的,有时候会不断的有新的测量数据加入插值节点集,
因此,通过寻找n个插值节点构造的的插值函数与n+1个插值节点构造的插值函数之间的关系,形成了牛顿插值法。推演牛顿插值法的方式是归纳法,也就是计算Ln(x)- Ln+1(x),并且从n=1开始不断的迭代来计算n+1时的插值函数。
牛顿插值法的公式是:
注意:在程序中我用W 代替
计算牛顿插值函数关键是要计算差商,n阶差商的表示方式如下:
关于差商我在这里并不讨论
计算n阶差商的公式是这样:
很明显计算n阶差商需要利用到两个n-1阶差商,这样在编程的时候很容易想到利用递归来实现计算n阶差商,不过需要注意的是递归有栈溢出的潜在危险,在计算差商的时候
更是如此,每一层递归都会包含两个递归,递归的总次数呈满二叉树分布:
这意味着递归次数会急剧增加:(。所以在具体的应用中还需要根据应用来改变思路或者优化代码。
废话少说,放码过来。
首先写最关键的一步,也就是计算n阶差商:
1 2 3 4 5 6 7 8 9 10 11 12 13 | """ @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 ]) |
看上面的牛顿插值函数公式,有了差商,还差
这个就比较好实现了:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | """ @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的闭包直接返回一个牛顿插值函数,闭包可以利用到它所处的函数之中的上下文数据。
1 2 3 4 5 6 7 8 9 10 11 | """ @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,这个从函数的表达式可以找到原因。
构造一些插值节点
1 2 3 | ''' 插值节点, 这里用二次函数生成插值节点,每两个节点x轴距离位10 ''' sr_x = [i for i in range ( - 50 , 51 , 10 )] sr_fx = [i * * 2 for i in sr_x] |
获得牛顿插值函数
1 2 3 4 | 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] # 根据插值函数获得测试用例的纵坐标 |
完整代码:
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 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 | # -*- 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() |
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步