import math x0 def lambda 100 y-x 1-x Rozenbrok lambda if not and not

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
import math
x0 = (2,2)
def F(x='',y=''):
k = 4
#f = lambda x,y: 100*(y-x**2)**2 + (1-x)**2 # Rozenbrok
f = lambda x,y: x**4 + 2*k*x**3 - k**2*x**2 - 2*k**3*x + y**4 - y**2*(2*k**2 + 2*k) + 2*k**4 + k**3 + k**2
if not x == '' and not y == '':
return f(x,y)
elif x == '' and not y == '':
return lambda x: f(x,y)
elif y == '' and not x == '':
return lambda y: f(x,y)
else:
return lambda x,y: f(x,y)
def gF(X,Y):
K = 4
#return [(F(x+1,y)-F(x,y))/1,(F(x,y+1)-F(x,y))/1]
return [4*X**3 + 2*K*X*X - 2*K*K*X - 2*K**3,4*Y**3 - 4*K*K*Y - 4*K*Y]
#return [400*X*Y+400*X**3+2*X-2,200*Y-2*X*X]
def min_func(f,point,step):
while True:
if f(point) < f(point+step):
break
point += step
return point
def search_min_of_func(f,point,step):
if f(point+step) < f(point-step):
return min_func(f,point,step)
elif f(point+step) == f(point-step):
if f(point+step) > f(point):
return point
else:
return min_func(f,point,step)
else:
return min_func(f,point,step*(-1))
x1 = [x0[0],x0[1]]
step = 0.0001
while True:
y = search_min_of_func(F(x=x1[0]),x1[1],step)
x = search_min_of_func(F(y=y),x1[0],step)
print x,y
if math.sqrt((x1[0]-x)**2+(x1[1]-y)**2) < 0.00000001:
print x,y
print F(x,y)
break
x1[0] = x
x1[1] = y
def gs(F,X,Y,Alfa,Step,Alfa_first):
[start_x,start_y] = [X,Y]
while True:
if F(X,Y) > F(X+Step*math.cos(Alfa),Y+Step*math.sin(Alfa)):
X = X+Step*math.cos(Alfa)
Y = Y+Step*math.sin(Alfa)
else:
break
if [start_x,start_y] == [X,Y]:
if abs(Alfa)-math.pi > abs(Alfa_first): return [X,Y]
else: return gs(F,X,Y,Alfa+math.pi/12,Step,Alfa_first)
else:
return [X,Y]
def get_alfa(X,Y):
if (X > 0 and Y > 0) or (X > 0 and Y < 0): return math.atan(Y/X)
else: return math.atan(Y/X)+math.pi
def result_gs(X,Y):
while True:
[Gx,Gy] = map(lambda A: -1*A,gF(X,Y))
Alfa = get_alfa(Gx,Gy)
[New_x,New_y] = gs(F(),X,Y,Alfa,0.0001,Alfa)
[New_gx,New_gy] = gF(New_x,New_y)
Grad = math.sqrt((New_x-X)**2+(New_y-Y)**2)
if Grad < 0.000001:
break
else:
[X,Y] = [New_x,New_y]
print [X,Y],F(X,Y)
print X,Y
result_gs(1,1)