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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
|
comment
Minimization
endcomment
function fmin (fff,a,b)
## Compute the Minimum of the convex function fff in [a,b],
## using the golden cut method.
## Additional parameters are passed to fff.
## fff may be an expression in x.
## You can specify an epsilon eps with eps=... as last parameter.
if (isvar("eps")); localepsilon(eps); endif;
x0=a; x3=b;
if isfunction(fff)
y0=fff(x0,args()); y3=fff(x3,args());
else
y0=expreval(fff,x0); y3=expreval(fff,x3);
endif
l=(3-sqrt(5))/2;
x1=x0+l*(x3-x0); x2=x3-l*(x3-x0);
if isfunction(fff)
y1=fff(x1,args()); y2=fff(x2,args());
else
y1=expreval(fff,x1); y2=expreval(fff,x2);
endif
repeat
if y1>y2;
x0=x1; x1=x2; x2=x3-l*(x3-x0);
y0=y1; y1=y2;
if isfunction(fff); y2=fff(x2,args());
else y2=expreval(fff,x2);
endif;
else
x3=x2; x2=x1; x1=x0+l*(x3-x0);
y3=y2; y2=y1;
if isfunction(fff); y1=fff(x1,args());
else y1=expreval(fff,x1);
endif;
endif;
if x0~=x3; return x0; endif;
end;
endfunction
function fmax (fff,a,b)
## Compute the Maximum of the concave function fff in [a,b],
## using the golden cut method.
## Additional parameters are passed to fff.
## fff may be an expression in x.
## You can specify an epsilon eps with eps=... as last parameter.
if (isvar("eps")); localepsilon(eps); endif;
x0=a; x3=b;
if isfunction(fff)
y0=fff(x0,args()); y3=fff(x3,args());
else
y0=expreval(fff,x0); y3=expreval(fff,x3);
endif
l=(3-sqrt(5))/2;
x1=x0+l*(x3-x0); x2=x3-l*(x3-x0);
if isfunction(fff)
y1=fff(x1,args()); y2=fff(x2,args());
else
y1=expreval(fff,x1); y2=expreval(fff,x2);
endif
repeat
if y1<y2;
x0=x1; x1=x2; x2=x3-l*(x3-x0);
y0=y1; y1=y2;
if isfunction(fff); y2=fff(x2,args());
else y2=expreval(fff,x2);
endif;
else
x3=x2; x2=x1; x1=x0+l*(x3-x0);
y3=y2; y2=y1;
if isfunction(fff); y1=fff(x1,args());
else y1=expreval(fff,x1);
endif;
endif;
if x0~=x3; return x0; endif;
end;
endfunction
function fextrema (fff,a,b,n=100)
## Find all internal extrema of fff in [a,b].
## fff may be an expression in x or a function.
## You can specify an epsilon eps with eps=... as last parameter.
## Returns {minima,maxima} (maybe of lengths 0)
if (isvar("eps")); localepsilon(eps); endif;
x=linspace(a,b,n);
mi=[]; ma=[];
if isfunction(fff);
y=fff(x;args());
else
y=expreval(fff,x);
endif;
loop 2 to n;
if y[#]>=y[#-1] && y[#]>y[#+1];
ma=ma|fmax(fff,x[#-1],x[#+1];args());
elseif y[#]<=y[#-1] && y[#]<y[#+1]
mi=mi|fmin(fff,x[#-1],x[#+1];args());
endif;
end;
return {mi,ma}
endfunction
function fffexprv (x,expr)
return evaluate(expr);
endfunction
function brentmin (fff,a,d=0.1,eps=epsilon())
## Compute the minimum of fff around d with Brent's method.
## d is an initial step size to look for a starting interval.
## eps is the final accuracy.
## Returns the point of minimal value.
## fff may be an expression in x.
## Return the point of minimum.
if isfunction(fff);
return brent(fff,a,d,eps;args());
else;
return brent("fffexprv",a,d,eps,fff);
endif;
endfunction
function neldermin (fff,v,d=0.1,eps=epsilon())
## Compute the minimum of fff around d with Nelder-Means's method.
## d is an initial step size for a starting simplex.
## eps is the final accuracy.
## Returns the point of minimal value.
## x is a 1xn vector.
## fff may be an expression in x (n>2).
## Return the point of minimum (1xn vector).
if isfunction(fff);
return nelder(fff,v,d,eps,args());
else;
return nelder("fffexprv",v,d,eps,fff);
endif;
endfunction
|
