/
Needle.m
143 lines (123 loc) · 4.75 KB
/
Needle.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
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
138
139
140
141
142
143
%% Needle on a corrugated surface
% Nick Trefethen and Hrothgar, December 2013
%%
% (Chebfun example opt/Needle.m)
% [Tags: #optimization]
function needle()
%%
% The final problem for Oxford's Numerical Analysis Problem Solving Squad this
% year was the following:
% <blockquote>
% A needle of length $1$ rests on a surface defined by the height function
% $$ h(x) = 0.1 x^2 + 0.1\sin(6x) + 0.03\sin(12x). $$
% What is the lowest possible height of the centre of the needle?
% </blockquote>
%%
% Let's begin with a picture:
LW = 'linewidth'; FS = 'fontsize'; MS = 'markersize';
s = chebfun('s',[-4 4]);
h = .1*s.^2 + .1*sin(6*s)+.03*sin(12*s);
close all, plot(h,LW,1), axis equal, axis([-4 4 -.4 2])
%%
% We see immediately that the optimal position of
% the center of the needle, call it $x$, will lie
% in $[-2,2]$. Actually it's pretty clear it will
% lie in $[-1,1]$.
%%
% This is obviously an optimization problem, but exactly how should we
% formulate it? One interpretation is that it is a problem of _semiinfinite
% programming_, because it mixes an objective function to minimize with a
% continuum of constraints. Presumably there are methods that could be used
% to solve it in this framework.
%%
% Here in the Oxford Numerical Analysis Group, when
% we hear the word "continuum", we think Chebfun. Suppose we
% specify two variables: $x$, the horizontal position of the center of the
% needle, and $\theta$, its angle counterclockwise from the horizontal. Given
% $x$ and $\theta$, we then ask how low the needle can lie. Let
%%
% $$ y(x,\theta) $$
%%
% be its minimal height, given that it does not cut below the surface. Then
% $y(x,\theta)$ is just a maximum of a continuous function over an interval,
% which we can compute with Chebfun like this:
function y = minfun(x,theta)
r = .5*cos(theta); hx = h{x-r,x+r};
needle = chebfun(@(s) tan(theta)*(s-x),[x-r x+r]);
y = max(hx - needle);
end
%%
% And here is a function for plotting a particular configuration:
function plotneedle(x,theta)
y = minfun(x,theta);
r = .5*cos(theta); hx = h{x-r,x+r};
needle = chebfun(@(s) y+tan(theta)*(s-x),[x-r x+r]);
hold off, plot(h,'b',needle,'k',LW,1)
axis equal, axis([-4 4 -.4 2])
end
%%
% For example, here are needle positions for $(x,\theta) = (-.6, -.2)$ and
% $(x,\theta) = (1.7, 1)$.
subplot(2,1,1)
plotneedle(-0.6,-0.2), title('needle with (x,theta) = (-0.6, -0.2)',FS,14)
subplot(2,1,2)
plotneedle(1.7,1), title('needle with (x,theta) = (1.7, 1)',FS,14)
%%
% Now we just have to minimize over $x$ and $\theta$.
% Let's first do that over a wide range.
npts = 25;
tic, x = linspace(-2,2,npts); theta = linspace(-1.5,1.5,npts);
[xx,thth] = meshgrid(x,theta); yy = 0*xx;
for k = 1:length(x)
for j = 1:length(theta)
yy(j,k) = minfun(xx(j,k), thth(j,k));
end
end
xxp = linspace(-2,2,100); ttp = linspace(-1.5,1.5,100)';
yyp = interp2(xx,thth,yy,xxp,ttp,'cubic');
close, contour(xxp,ttp,yyp,80), grid on, xlabel('x',FS,14), ylabel('theta',FS,14)
colorbar, title(['min value on grid: ' num2str(min(yy(:)))],FS,14), toc
%%
% In this picture we see that there are two promising
% regions: one with $(x,\theta) \approx (-.5, -.4)$,
% and one with $(x,\theta) \approx (.5, -.2)$.
% The central white regions have an interesting interpretation: if the needle
% is balanced on top of a mountain, then moving it left or right, or tilting
% it, doesn't have much effect.
%%
% Zooming in confirms this picture:
tic, x = linspace(-0.8,0.6,npts); theta = linspace(-0.5,0,npts);
[xx,thth] = meshgrid(x,theta); yy = 0*xx;
for k = 1:length(x)
for j = 1:length(theta)
yy(j,k) = minfun(xx(j,k), thth(j,k));
end
end
xxp = linspace(-0.8,0.6,100); ttp = linspace(-0.5,0,100)';
yyp = interp2(xx,thth,yy,xxp,ttp,'cubic');
levels = 0.06:.003:0.12;
close, contour(xxp,ttp,yyp,levels), grid on, xlabel('x',FS,14), ylabel('theta',FS,14)
title(['min value on grid: ' num2str(min(min(yy)))],FS,14), toc
%%
% The winner seems to be the region on the right. From here the right thing to
% do is call a bivariate optimization routine. In basic MATLAB the simplest
% one is the direct search code
% <a href='http://www.mathworks.co.uk/help/matlab/ref/fminsearch.html'>fminsearch</a>.
% This requires the input to be a single vector, so we'll need a wrapper:
function y = minfunwrapper(xvec)
y = minfun(xvec(1), xvec(2));
end
%%
% Here goes.
opts = optimset('tolx',1e-14,'display','off');
guess = [.41, -0.2];
tic, [xvec,yval] = fminsearch(@minfunwrapper,guess,opts); toc
%%
% So it would seem that to 10 digits or more, the minimal height is around
yval
%%
% Here is a closeup of the solution:
plotneedle(xvec(1),xvec(2)), hold on
axis equal, axis([-2 2 -.4 1.2])
plot(xvec(1),yval,'.k',MS,12), grid on
end