Третий алгоритм Пирсона для поиска минимума

import math
import numpy as np
import sympy as sym
 
delta = 1E-3
eps_find = 1E-7
 
def f_Square(x1, x2):
    return 100*(x2-x1)**2 + (1-x1)**2
 
def df_Square(X):
    x1, x2 = sym.symbols('x1 x2')
    c1 = (sym.diff(f_Square(x1, x2),x1)).subs(x1,X.item(0)).subs(x2,X.item(1))
    c2 = (sym.diff(f_Square(x1, x2),x2)).subs(x1,X.item(0)).subs(x2,X.item(1))
    return np.matrix(str(c1)+";"+str(c2))
 
def ddf_Square(X):
    x1, x2 = sym.symbols('x1 x2')
    c1 = sym.diff(f_Square(x1, x2),x1)
    c2 = sym.diff(f_Square(x1, x2),x2)
    c11 = (sym.diff(c1,x1)).subs(x1,X.item(0)).subs(x2,X.item(1))
    c12 = (sym.diff(c1,x2)).subs(x1,X.item(0)).subs(x2,X.item(1))
    c21 = (sym.diff(c2,x1)).subs(x1,X.item(0)).subs(x2,X.item(1))
    c22 = (sym.diff(c2,x2)).subs(x1,X.item(0)).subs(x2,X.item(1))
    return np.matrix(str(c11)+","+str(c12)+";"+str(c21)+","+str(c22))
 
def f_Rose(x1, x2):
    return 100*(x2-x1**2)**2 + (1-x1)**2
 
def df_Rose(X):
    x1, x2 = sym.symbols('x1 x2')
    c1 = sym.diff(f_Rose(x1, x2),x1).subs(x1,X.item(0)).subs(x2,X.item(1))
    c2 = sym.diff(f_Rose(x1, x2),x2).subs(x1,X.item(0)).subs(x2,X.item(1))
    return np.matrix(str(c1)+";"+str(c2))
 
def ddf_Rose(X):
    x1, x2 = sym.symbols('x1 x2')
    c1 = sym.diff(f_Rose(x1, x2),x1)
    c2 = sym.diff(f_Rose(x1, x2),x2)
    c11 = (sym.diff(c1,x1)).subs(x1,X.item(0)).subs(x2,X.item(1))
    c12 = (sym.diff(c1,x2)).subs(x1,X.item(0)).subs(x2,X.item(1))
    c21 = (sym.diff(c2,x1)).subs(x1,X.item(0)).subs(x2,X.item(1))
    c22 = (sym.diff(c2,x2)).subs(x1,X.item(0)).subs(x2,X.item(1))
    return np.matrix(str(c11)+","+str(c12)+";"+str(c21)+","+str(c22))
 
def f_Test(x, y):
    return -(2*math.exp(-((x-1)/2)**2-((y-1)/2)**2)+3*math.exp(((x-1)/3)**2+((y-3)/3)**2))
 
def df_Test(X):
    x, y = sym.symbols('x y')
    c1 = (sym.diff(f_Test(x, y),x)).subs(x,X.item(0)).subs(y,X.item(1))
    c2 = (sym.diff(f_Test(x, y),y)).subs(x,X.item(0)).subs(y,X.item(1))
    return np.matrix(str(c1)+";"+str(c2))
 
def ddf_Test(X):
    x, y = sym.symbols('x1 x2')
    c1 = sym.diff(f_Test(x, y),x)
    c2 = sym.diff(f_Test(x, y),y)
    c11 = (sym.diff(c1,x)).subs(x,X.item(0)).subs(y,X.item(1))
    c12 = (sym.diff(c1,y)).subs(x,X.item(0)).subs(y,X.item(1))
    c21 = (sym.diff(c2,x)).subs(x,X.item(0)).subs(y,X.item(1))
    c22 = (sym.diff(c2,y)).subs(x,X.item(0)).subs(y,X.item(1))
    return np.matrix(str(c11)+","+str(c12)+";"+str(c21)+","+str(c22))
 
def mini(f, x_prev):
    d = delta
 
    a, b, x_curr, x_next = 0, 0, 0, 0
 
    if (f(x_prev) > f(x_prev+d)):
        x_curr = x_prev + d
        h = d
    elif (f(x_prev) > f(x_prev - d)):
        x_curr = x_prev - d
        h = -delta
    else:
        return x_prev - d, x_prev + d
 
    while (True):
        h *= 2
        x_next = x_curr + h
 
        if (f(x_curr) > f(x_next)):
            x_prev = x_curr
            x_curr = x_next
        else:
            a = x_prev
            b = x_next
            if (a > b):
                a, b = b, a
            return a, b
 
def getFn(n):
    return (int)((((1+5**0.5)/2)**n - ((1-5**0.5)/2)**n)/5**0.5)
 
def fibo(f, e, a, b, out = False):
    i = 1
 
    n = 1
 
    if (b < a):
        a, b = b, a
 
    while (getFn(n+2)<(b-a)/e):
        n += 1
 
    n-=1
 
    ap = a
    bp = b
 
    x1 = a + getFn(n)*(b - a) / getFn(n+2)
    x2 = a + getFn(n+1)*(b - a) / getFn(n+2)
 
    fx1 = f(x1)
    fx2 = f(x2)
 
    while (True):
        if (i > n+1):
            break
 
        ap = a
        bp = b
 
        if (fx1 < fx2):
            b = x2
            x2 = x1
            fx2 = fx1
            x1 = a + getFn(n-i+1)*(b - a) / getFn(n-i+3)
            fx1 = f(x1)
        else:
            a = x1
            x1 = x2
            fx1 = fx2
            x2 = a + getFn(n-i+2)*(b - a) / getFn(n-i+3)
            fx2 = f(x2)
 
        i += 1
    return (a+b)/2
 
def parabol(f, a, b):
    x = (a + b) / 2
    h = delta
    fx = f(x)
    fxh = f(x + h)
    fhx = f(x - h)
    i = 0
    tau = 0.5
 
    while (h > 2 * eps_find and (fxh - 2 * fx + fhx)!=0):
        xk = x - h / 2 * (fxh - fhx) / (fxh - 2 * fx + fhx)
        if (fxh - 2 * fx + fhx <= 0):
            h = -3 * h
        else:
            fxk = f(xk)
            if (fxk >= fx):
                h = tau * (xk - x)
            else:
                h = abs(x - xk) / 2
                x = xk
                fx = fxk
                fxh = f(x + h)
                fhx = f(x - h)
        i += 1
 
    return x
 
### Метод Пирсона
n=[]
x_=[]
dx=[]
dg=[]
lam=[]
n.append(np.identity(2))
x_.append(np.matrix("-3;3"))
 
funcs = [[f_Square,df_Square,ddf_Square],
         [f_Rose,df_Rose,ddf_Rose],
         [f_Test,df_Test,ddf_Test]]
 
mfun = funcs[1]
for i in range(1,10):
    k = i-1
 
    flam = lambda lem: mfun[0](
        (x_[k]-lem*n[k]*mfun[1](x_[k]))[0],
        (x_[k]-lem*n[k]*mfun[1](x_[k]))[1]
    ).item(0)
 
    #print(str(x_[k])+"-lem*"+str(n[k])+"*"+str(mfun[1](x_[k]))+"\n")
 
    argminflam = mini(flam, 0)
 
    #fminimumlam = fibo(flam, 0.000001, argminflam[0], argminflam[1])
    fminimumlam = parabol(flam, argminflam[0], argminflam[1])
 
    lam.append(fminimumlam)
 
    print(fminimumlam, flam(fminimumlam), "\n")
 
    x_.append(x_[k]-lam[k]*n[k]*mfun[1](x_[k]))
 
    print(argminflam[0], argminflam[1], "\n")
    #print(mfun[0](x_[i][0], x_[i][1]), x_[i])
    print(mfun[0](x_[k][0], x_[k][1]), x_[k])
 
    dx.append(x_[i]-x_[k])
    dg.append(mfun[1](x_[i])-mfun[1](x_[k]))
    t1 = (dx[k]-n[k]*dg[k])
    t2 = (n[k]*dg[k]).T
    t3 = (dg[k].T*n[k]*dg[k]).item(0)
 
    print("========")
    if (t3 == 0):
        break
 
    n.append(n[k] + (t1*t2)/t3)
k+=1
 
print(x_[k])