import math


def AdamusExplicitly(k,first_ys,x0,y0,h,xk,fxy,fx_real):
    betas = [
        [1],
        [3/2, -1/2],
        [23/12, -16/12, 5/12],
        [55/24, -59/24, 37/24, -9/24],
        [1901/720, -2774/720, 2626/720, -1274/720, 251/720],
        [4277/1440, -7923/1600, 9982/1440, -7298/1440, 2877/1440, -476/1440]
    ]
    # Bks = [1/2,5/12,3/8,251/720,95/288,10987/60480]
    result = []
    if k<0 or k >= len(betas):
        print("k超出范围")
        return None
    if len(first_ys) != k + 1:
        print("first_ys(前几个y的值)与k长度不匹配")
        return None
    fxys = [fxy(x0+i*h, first_ys[i]) for i in range(k + 1)]
    for i in range(k + 1):
        result.append((x0 + i * h, first_ys[i]))
    x1 = x0 + h*k
    y0 = first_ys[k]
    while x1 < xk:
        x1 += h
        delta_y = h * sum(betas[k][i] * fxys[k-i] for i in range(k + 1))
        y1 = y0 + delta_y
        print(f"x={x1}, y={y1}, delta_y={delta_y}, 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
    return result



if __name__ == "__main__":
    # 定义初始值和参数#####################################################################################
    k = 3                    # 阿达姆斯k+1步显式方法，精度为k+1阶，最常用k=3，其他阶数我没试过   P253
    first_ys = [1, 0.904837418036, 0.818730753078, 0.7408182206817]  # 前几个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 = AdamusExplicitly(k, first_ys, x0, y0, h, xk, fxy, fx_real)
    print("计算结果看到有几.99999或者几.00000就自己四舍五入一下，有可能会多算一点，自己比较一下")
    if result:
        for xy in result:
            print(f"(x, y): {xy}")