-
Notifications
You must be signed in to change notification settings - Fork 0
/
weiguo_new.m
107 lines (103 loc) · 2.04 KB
/
weiguo_new.m
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
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
% 黄金分割法
% initial environment
a1 = -0.01;
a4 = 0.01;
epsilon = 0.0001;
a2 = a1 + 0.382*(a4 - a1);
x = [3 7]';
S = [0.7 0.714]';
x2 = x + a2 * S;
temp = x;
x = x2;
f2 = x(1)^5 + x(2)^4 -5*x(1)^4+30*x(2)^3 +25*x(1)^3+10*x(2)^2+x(1)*x(2) - 10*x(1) - 14*x(2) + 30;
x = temp;
a3 = a1 + 0.618*(a4 - a1);
x3 = x + a3 * S;
x = x3;
f3 = x(1)^5 + x(2)^4 -5*x(1)^4+30*x(2)^3 +25*x(1)^3+10*x(2)^2+x(1)*x(2) - 10*x(1) - 14*x(2) + 30;
x = temp;
% global varies
astar = 0;x1 = 0;fstar = 0;
while 1
if f2 < f3
a4 = a3;a3 = a2;f3 = f2;
a2 = a1 + 0.382*(a4 - a1);
x2 = x + a2*S;
x = x2;
f2 = x(1)^5 + x(2)^4 -5*x(1)^4+30*x(2)^3 +25*x(1)^3+10*x(2)^2+x(1)*x(2) - 10*x(1) - 14*x(2) + 30;
x = temp;
else
a1 = a2;a2 = a3;f2 = f3;
a3 = a1 + 0.618*(a4 - a1);
x3 = x + a3 * S;
x = x3;
f3 = x(1)^5 + x(2)^4 -5*x(1)^4+30*x(2)^3 +25*x(1)^3+10*x(2)^2+x(1)*x(2) - 10*x(1) - 14*x(2) + 30;
x = temp;
end
if abs(a4 - a1) > epsilon
continue;
else
astar = (a1 + a4)/2;
a1,a4
x1 = x + astar * S;
x = x1;
fstar = x(1)^5 + x(2)^4 -5*x(1)^4+30*x(2)^3 +25*x(1)^3+10*x(2)^2+x(1)*x(2) - 10*x(1) - 14*x(2) + 30;
x = temp;
break;
end
end
astar,fstar
% 伟哥自己的二次插值修改
e=0.001;
while 1
s1 = (f3 - f2) / (a3 - a2);
s2 = ((f2 - f1) / (a2 - a1) - s1) / (a2 - a3);
if s2==0
f1 = f2;
x1 = x + a2 * s;
break
else
a0 = (a1 + a3 - s1 / s2) / 2;
if (a0 - a1) * (a3 - a0) < 0
f1 = f2;
x1 = x + a2 * s;
break
else
x0 = x + a4 * s;
f0 = x0(1)^5 + x0(2)^4 -5 * x0(1)^4+30*x0(2)^3 +25*x0(1)^3+10*x0(2)^2+x0(1)*x0(2) - 10*x0(1) - 14*x0(2) + 30;
if abs(a0 - a2) < e
if f0<f2
f1 = f0;
x1 = x0;
break
else
f1 = f2;
x1 = x + a2 * s;
break
end
else
if (a0-a2) < 0
if f2 < f0
a1 = a0;
f1 = f0;
continue
else
a3 = a2; f3 = f2;
a2 = a0; f2 = f0;
continue
end
else
if f2 < f0
a3 = a0;
f3 = f0;
continue
else
a1 = a2; f1 = f2;
a2 = a0; f2 = f0;
continue
end
end
end
end
end
end