import math


def AdamusImplicitly(k,first_ys,x0,y0,h,xk,fxy,fx_real):
    betas = [
        [1],
        [1/2, 1/2],
        [5/12, 8/12, -1/12],
        [9/24, 19/24, -5/24, 1/24],
        [251/720, 646/720, -261/720, 106/720, -19/720],
        [475/1440, 1427/1440, -798/1440, 482/1440, -173/1440, 27/1440]
    ]
    
    result = []
    if k<0 or k >= len(betas):
        print("k超出范围")
        return None
    if len(first_ys) != k:
        print("first_ys(前几个y的值)与k长度不匹配")
        return None
    fxys = [fxy(x0+i*h, first_ys[i]) for i in range(k)]
    for i in range(k):
        result.append((x0 + i * h, first_ys[i]))
    x1 = x0 + h*k
    y0 = first_ys[k-1]
    while x1 < xk:
        y1_t = y0
        y1 = y0 + 1
        for j in range(1000):
            y1 = y0 + h * sum(betas[k][i+1] * fxys[k-i-1] for i in range(k))+ h * betas[k][0] * fxy(x1, y1_t)
            if abs(y1 - y1_t) < 1e-14:
                break
            y1_t = y1

        # y1 = y0 + delta_y
        print(f"x={x1}, y={y1}, real_y={fx_real(x1)}, abs(real_y-y)={abs(fx_real(x1) - y1)}")
        result.append((x1, y1))
        fxys.append(fxy(x1, y1))
        fxys.pop(0)
        y0 = y1
        x1 += h
    return result



if __name__ == "__main__":
    # 定义初始值和参数#####################################################################################
    k = 3                    # 阿达姆斯k步隐式方法，精度为k+1阶，最常用k=3，其他阶数我没试过   P255
    first_ys = [1, 0.904837418036, 0.818730753078]  # 前几个y的值，可用龙格-库塔计算或者知道精确解自己算出来
    x0 = 0.0                 # 初始x值
    y0 = 1.0                 # 初始y值
    h = 0.1                  # 步长
    xk = 1.0                 # 最终x值
    fxy = lambda x, y: -y                   # 定义f(x,y)函数，dx/dy = f(x,y)，导函数
    fx_real = lambda x: math.exp(-x)        # 实际解函数，用于验证结果，如果不知道或者不用算误差，可以直接写个 lambda x: 0
    # 调用阿达姆斯隐式方法
    result = AdamusImplicitly(k, first_ys, x0, y0, h, xk, fxy, fx_real)
    print("计算结果看到有几.99999或者几.00000就自己四舍五入一下，有可能会多算一点，自己比较一下")
    if result:
        for xy in result:
            print(f"(x, y): {xy}")