''' mbinary
#########################################################################
# File : numerical_integration.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 : numerical integration.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.xyz
# Github: https://github.com/mbinary
# Created Time: 2018-05-11  08:58
# Description:
#       numerical intergration: using Newton-Cotes integration,  and Simpson
#       数值积分, 使用 牛顿-科特斯积分, 辛普森
#########################################################################


import numpy as np


def trapezoidal(a, b, h, fs):
    '''梯形积分公式'''
    xs = [i for i in np.arange(a, b+h, h)]
    print(xs)
    ret = h*(sum(fs)-fs[0]/2 - fs[-1]/2)
    print(ret)
    return ret


def simpson(a, b, h, fs):
    '''辛普森积分公式'''
    xs = [i for i in np.arange(a, b+h, h)]
    print(xs)
    ret = h/3*(4 * sum(fs[1::2]) + 2*sum(fs[2:-1:2]) + fs[0]+fs[-1])
    print(ret)
    return ret


def romberg(a, b, f, epcilon):
    '''romberg(龙贝格) 数值积分'''
    h = b-a
    lst1 = [h*(f(a)+f(b))/2]
    print(lst1)
    delta = epcilon
    k = 1
    while delta >= epcilon:
        h /= 2
        k += 1
        lst2 = []
        lst2.append((lst1[0]+h*2*sum(f(a+(2*i-1)*h)
                                     for i in range(1, 2**(k-2)+1)))/2)
        for j in range(0, k-1):
            lst2.append(lst2[j]+(lst2[j]-lst1[j])/(4**(j+1)-1))
        delta = abs(lst2[-1]-lst1[-1])
        lst1 = lst2
        print(lst1)


if __name__ == '__main__':
    a, b, h = 0.6, 1.8, 0.2
    fs = [5.7, 4.6, 3.5, 3.7, 4.9, 5.2, 5.5]
    trapezoidal(a, b, h, fs)
    simpson(a, b, h, fs)
    romberg(1, 2, lambda x: sin(x**4), 1e-4)