【python代码】坐标轮换法

这里的坐标轮换法只更新二维的坐标，为了明确更新方法，这里没有用矩阵和向量，具体算法见相关资料或者下面代码的注释。（需要注意的是，我这里写的代码是寻找函数的极小值）

import math


# 定义极小值
eps = 1e-6


# ################## 坐标轮换法

# 函数方程 f = a * x_1^2 + b * x_2^2 + c * x_1 * x_2 + d * x_1 + e * x_2 + f
def get_func2(a, b, c, d, e, f):
    def func(x_1=None, x_2=None):
        return a * x_1 ** 2 + b * x_2 ** 2 + c * x_1 * x_2 + d * x_1 + e * x_2 + f
    return func


# 黄金分割法用于一维坐标的更新
def search_1d(func, input_range, x_1, x_2, dim=0):  # dim=0代表更新x_1,dim=1代表更新x_2，下面涉及到的dim也同理
    # 接受每次变化的区间
    intervals = [input_range]
    # 计算每次的区间大小，进行迭代
    dif = intervals[-1][1] - intervals[-1][0]
    while dif > eps:
        left = intervals[-1][1] - dif * 0.618
        right = intervals[-1][0] + dif * 0.618
        if dim == 0:
            f_l = func(x_1=left, x_2=x_2)
            f_r = func(x_1=right, x_2=x_2)
        else:
            f_l = func(x_1=x_1, x_2=left)
            f_r = func(x_1=x_1, x_2=right)
        if f_l < f_r:
            left = intervals[-1][0]
        elif f_l > f_r:
            right = intervals[-1][1]
        intervals.append([left, right])
        dif = right - left
    x_star = (intervals[-1][0] + intervals[-1][1]) / 2
    return x_star


# 进退法确定更新区间大小（进退法的原理可以参照罗尔定理，具体自行查阅相关博客）
def cal_intervals(func, x_1, x_2, h=0.1, dim=0):
    if dim == 0:
        if func(x_1 + h, x_2) - func(x_1, x_2) < eps:
            return x_1, x_1 + h

        if func(x_1 + h, x_2) > func(x_1, x_2):
            l = 1
            while func(x_1 - l * h, x_2) < func(x_1, x_2):
                l += 1
            return x_1 - l * h, x_1 + h

        l = 1
        r = 2
        while func(x_1 + r * h, x_2) < func(x_1, x_2):
            r += 1
        while func(x_1 - l * h , x_2) < func(x_1, x_2):
            l += 1
        return x_1 - l * h, x_1 + r * h

    else:
        if func(x_1, x_2 + h) - func(x_1, x_2) < eps:
            return x_2, x_2 + h

        if func(x_1, x_2 + h) > func(x_1, x_2):
            l = 1
            while func(x_1, x_2 - l * h) < func(x_1, x_2):
                l += 1
            return x_2 - l * h, x_2 + h

        l = 1
        r = 2
        while func(x_1, x_2 + r * h) < func(x_1, x_2):
            r += 1
        while func(x_1, x_2 - l * h) < func(x_1, x_2):
            l += 1
        return x_2 - l * h, x_2 + r * h


# 坐标轮换法
def coordinate_asc(func, x_1, x_2):
    # 记录上一次的坐标值
    last_x_1 = x_1
    last_x_2 = x_2
    while True:
        # 更新x_1
        x_1_l, x_1_r = cal_intervals(func, x_1, x_2, dim=0)  # 用进退法确定x_1更新区间
        x_1 = search_1d(func, [x_1_l, x_1_r], x_1, x_2, dim=0)  # 用黄金分割法更新x_1的值
        # 更新x_2
        x_2_l, x_2_r = cal_intervals(func, x_1, x_2, dim=1)  # 用进退法确定x_2更新区间
        x_2 = search_1d(func, [x_2_l, x_2_r], x_1, x_2, dim=1)  # 用黄金分割法更新x_2的值
        # 判断是否停止（如果更新后的坐标与上一次的坐标的距离小于等于极小值）
        if math.sqrt((x_1 - last_x_1) ** 2 + (x_2 - last_x_2) ** 2) <= eps:
            break
        # 记录上一次的坐标值
        last_x_1 = x_1
        last_x_2 = x_2
    return x_1, x_2


if __name__ == '__main__':
    # 输入确定方程的参数
    par = input("请输入参数a b c d e f: ")
    str = par.split(',')
    par_int = [int(i) for i in str]
    # 得到函数方程
    func2 = get_func2(par_int[0], par_int[1], par_int[2], par_int[3], par_int[4], par_int[5])
    # 坐标轮换法
    x_1_star, x_2_star = coordinate_asc(func2, 3, 3)