采用三次插值法进行一维搜索

简介

三次插值法是用a,b两点处的函数值$f(a),f(b)$和导数值$f^{{}^{\prime }}(a),f^{{}^{\prime }}(b)$来构造三次插值多项式，并以$g(x)$的极小点作为$f(x)$的近似极小点。为了保证$f(x)$的极小点$x^{\ast }\in [a,b]$，假定$f(x)\in C^1$，并满足$$a<b,f^{{}^{\prime }}(a)<0,f^{{}^{\prime }}(b)>0$$插值多项式为:$$g(x)=\alpha (x-a{)}^3+\beta (x-a{)}^2+\gamma (x-a)+\delta$$由插值条件可知$\alpha ,\beta ,\gamma ,\delta$应该满足下式：$$\{\begin{array}{rl}& \delta =f(a)\\ & \gamma =f^{{}^{\prime }}(a)\\ & \alpha (b-a{)}^3+\beta (b-a{)}^2+\gamma (b-a)+\delta =f(b)\\ & 3\alpha (b-a{)}^2+2\beta (b-a)+\gamma =f^{{}^{\prime }}(b)\end{array}$$
通过推导可得近似极小点$\overline{x}$的计算公式如下：$$\overline{x}=a-\frac{\gamma }{\beta +\sqrt{{\beta }^2-3\alpha \gamma }}$$
其中$\alpha ,\beta$的计算公式如下：$$\{\begin{array}{rl}& u=\frac{f(b)-\delta }{b-a}-\gamma \\ & v=f^{{}^{\prime }}(b)-\gamma \\ & \beta =\frac{3u-v}{b-a}\\ & \alpha =\frac{v-2u}{(b-a{)}^2}\end{array}$$
当$|f^{{}^{\prime }}(\overline{x})|<\epsilon$时，$\epsilon$为给定的收敛精度，$x^{\ast }\approx \overline{x}$；若$f^{{}^{\prime }}(\overline{x})>0$，用$\overline{x}$取代$b$；若$f^{{}^{\prime }}(\overline{x})<0$，用$\overline{x}$取代$a$,不断迭代直至满足收敛条件。

优化问题

同抛物线法相同，求$f(x)=x^4-4x^3-6x^2-16x+4$的极小值点。由于$f(0)=4,f^{{}^{\prime }}(0)=-16,f(10)=5244,f^{{}^{\prime }}(10)=2664$，因此，选取为$a=0,b=10$，收敛精度$\epsilon =1e-3$，计算程序如下：

import sys
def cal_fun_f(x):
    """
    函数值计算函数
    """
    return x**4-4*x**3-6*x**2-16*x+4

def derivative_cal_f(x):
    """
    一阶导数计算函数
    """
    return 4*x**3-12*x**2-12*x-16

def derivative_cal_f2(x):
    """
    二阶导数计算函数
    """
    return 12*x**2-24*x-12

def get_iteration_point(a, b):
    """
    得到迭代点
    :param a:
    :param b:
    """
    delta = cal_fun_f(a)
    gamma = derivative_cal_f(a)
    f_b = cal_fun_f(b)
    d_b = cal_fun_f(b)
    u = (f_b-delta)/(b-a)-gamma
    v = d_b-gamma
    beta = (3*u-v)/(b-a)
    alpha = (v-2*u)/(b-a)**2
    return a-gamma/(beta+(beta**2-3*alpha*gamma)**0.5)


def cubic_interpolation(a, b, gap):
    """
    采用三次插值法求解一维函数的极小值点
    :param gap: 收敛精度
    :param a: 区间左端点
    :param b: 区间右端点
    """
    # 检查是否满足假设
    if (a<b) and (derivative_cal_f(a) < 0) and (derivative_cal_f(b) > 0):
        pass
    else:
        sys.exit("输入条件不满足假设")
    i = 0
    x = get_iteration_point(a, b)
    print(f"第{i}次迭代极小值点", x)
    while abs(derivative_cal_f(x)) > gap:
        if derivative_cal_f(x) > 0:
            b = x
        else:
            a = x
        x = get_iteration_point(a, b)
        i += 1
        print(f"第{i}次迭代极小值点", x)

    # 检验二阶导数是否满足
    print(f"极小值点的二阶导数", derivative_cal_f2(x))

if __name__ == "__main__":
    cubic_interpolation(0, 10, 1e-3)

结果如下：

第31次迭代极小值点 3.9998833532632627
第32次迭代极小值点 3.999913181303018
第33次迭代极小值点 3.999935381702105
第34次迭代极小值点 3.999951905109375
第35次迭代极小值点 3.9999642032791956
第36次迭代极小值点 3.9999733566886855
第37次迭代极小值点 3.9999801695011494
第38次迭代极小值点 3.9999852402350826
第39次迭代极小值点 3.999989014356277
极小值点的二阶导数 83.99920903510015

可以看到成功找到了极小值点4.