algorithm-in-python/math/numericalAnalysis/interplotion.py

''' mbinary
#########################################################################
# File : interplotion.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.xyz
# Github: https://github.com/mbinary
# Created Time: 2018-10-02  21:14
# Description:
#########################################################################
'''

#########################################################################
# File : interplotion.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.xyz
# Github: https://github.com/mbinary
# Created Time: 2018-05-18  09:29
# Description: 插值计算,有牛顿插值,拉格朗日插值,以及通过插值得到的多项式估计新的函数值
#########################################################################


import sympy
from collections import namedtuple
from functools import reduce
from operator import mul

X = sympy.Symbol('x')
point = namedtuple('point', ['x', 'y'])


class interplotion:
    def __init__(self, points):
        self.points = [point(x, y) for x, y in points]
        self.xs = [i for i, j in points]
        self.poly, self.rem = self.newton(self.points, 0, len(self.points)-1)

    def newton(self, li, a, b):
        '''li:[(x,f(x))...]'''

        qs = [li[0].y]

        def quoDiff(begin, end):
            if begin == end:
                return li[begin].y
            q = (quoDiff(begin+1, end)-quoDiff(begin, end-1)) / \
                (li[end].x-li[begin].x)
            if begin == a:
                qs.append(q)
            return q

        quoDiff(a, b)
        poly, base = 0, 1
        for i, q in enumerate(qs):
            poly += q*base
            base *= X-li[i].x
        return poly, base*qs[-1]

    def lagrange(self, points=None):
        xs = None
        if points is None:
            xs = self.xs
            points = self.points
        else:
            xs = [x for x, y in points]
        product = reduce(mul, [X-x for x in xs], 1)
        poly = 0
        for x, y in points:
            tmp = product/(X-x)
            coef = y/(tmp.subs(X, x))
            poly += coef * tmp
        return poly

    def predict(self, val, poly=None):
        if poly is None:
            poly = self.poly
        return poly.subs(X, val)  # note the func  subs


if __name__ == '__main__':
    f = interplotion([(81, 9), (100, 10), (121, 11)])
    p = f.lagrange()
    print(p.subs(X, 105))
    print(p)

    intor = interplotion([(0, 11), (0.02, 9), (0.04, 7), (0.06, 10)])
    p = intor.lagrange()
    print(p)
    res = intor.predict(0.08)
    print(res)