colloc.py 文件源码

def calc_cp(alpha,beta,k):
    gamma = []
    for i in range(k+1):
        num = factorial(alpha+k)*factorial(alpha+beta+k+i)
        denom = factorial(alpha+i)*factorial(k-i)*factorial(i)
        gamma.insert(i,num/denom)

    poly = []
    for i in range(k+1):
        if i == 0:
            poly.insert(i,gamma[i])
        else:
            prod = [1]
            j=1
            while j<=i:
                prod=conv(prod,[1,-1])
                j=j+1
            while len(poly)<len(prod):
                poly.insert(0,0)
            prod = [gamma[i]*t for t in prod]
            poly = [sum(pair) for pair in zip(poly,prod)]

    cp = numpy.roots(poly)
    return cp