双约束重力模型

c = [[0, 0, 0, 0], [0, 2, 4, 4], [0, 4, 1, 2], [0, 4, 2, 2]]
o = [0, 8, 14, 10]
d = [0, 8, 14, 10]
am = [0, 1, 1, 1]
ampl1 = [0, 0, 0, 0]
bm = [0, 0, 0, 0]
bmpl1 = [0, 0, 0, 0]
mi = 0.01
e = 0.05
m = 0
for j in range(1, 4):
    sumnum = 0
    for i in range(1, 4):
        sumnum += am[i] * o[i] * pow(c[i][j], -mi)
    bm[j] = 1 / sumnum
for i in range(1, 4):
    sumnum = 0
    for j in range(1, 4):
        sumnum += bm[j] * d[j] * pow(c[i][j], -mi)
    ampl1[i] = 1 / sumnum
for j in range(1, 4):
    sumnum = 0
    for i in range(1, 4):
        sumnum += ampl1[i] * o[i] * pow(c[i][j], -mi)

    bmpl1[j] = 1 / sumnum
def ifif(a, b, c, d):
    f = [0, 0, 0, 0, 0]
    for i in range(1, 4):
        if (a[i] / b[i] > 1 - e and a[i] / b[i] < 1 - e) and (c[i] / d[i] > 1 - e and c[i] / d[i] < 1 + e):
            f[i] = True
    if (f[1] and f[2]) and f[3]:
        return True
    else:
        return False
for k in range(50):
    if not ifif(ampl1, am, bmpl1, bm):
        m += 1
        mi += 0.05
        for j in range(1, 4):
            sumnum = 0
            for i in range(1, 4):
                sumnum += am[i] * o[i] * pow(c[j][i], -mi)
            bm[j] = 1 / sumnum
            print(bm[j])
        for i in range(1, 4):
            sumnum = 0
            for j in range(1, 4):
                sumnum += bm[j] * d[j] * pow(c[i][j], -mi)
            ampl1[i] = 1 / sumnum
        for j in range(1, 4):
            sumnum = 0
            for i in range(1, 4):
                sumnum += ampl1[i] * o[i] * pow(c[j][i], -mi)
            bmpl1[j] = 1 / sumnum
print(am, ampl1, bm, bmpl1, mi)