Problem 4:

What is the global minimum of the function

exp(sin(50*x))+sin(60*exp(y))+sin(70*sin(x))+sin(sin(80*y))-sin(10*(x+y))+(x^2+y^2)/4 ?

Solution:

Look at the function for reasonably small values of x and y.

>    restart;
F := (x,y) -> exp(sin(50*x)) + sin(60*exp(y)) + sin(70*sin(x)) +
sin(sin(80*y)) - sin(10*(x+y)) + 1/4*(x^2+y^2);

F := proc (x, y) options operator, arrow; exp(sin(50*x))+sin(60*exp(y))+sin(70*sin(x))+sin(sin(80*y))-sin(10*x+10*y)+1/4*x^2+1/4*y^2 end proc
F := proc (x, y) options operator, arrow; exp(sin(50*x))+sin(60*exp(y))+sin(70*sin(x))+sin(sin(80*y))-sin(10*x+10*y)+1/4*x^2+1/4*y^2 end proc

>    plot3d(F(x,y),x=-2..2,y=-2..2,axes=box,grid=[100,100],style=patchcontour,shading=ZHUE);

[Maple Plot]

It looks like a nightmare for someone trying to find the global minimum: there are huge numbers of local minima.  Note however that there's an easy lower bound: exp(-1)-4+(x^2+y^2)/4 <= F(x,y) .  From the graph it looks like there are points where F(x,y) is below -2.  In fact, the minimum value plotted is

>    P:= %:

>    G:= op([1,3],P);

G := Array(%id = 28615388)

>    CurrMin:=min(seq(seq(G[i,j],i=1..100),j=1..100));

CurrMin := -2.81950371067638050

For any point doing better than that, exp(-1)-4+(x^2+y^2)/4 < CurrMin  so x^2+y^2 < 4*(CurrMin+4-exp(-1))  

>    rmax:=sqrt(4*(CurrMin+4-exp(-1.0)));

rmax := 1.802905265

So we certainly don't have to search farther away from the origin than that.  Now suppose we search in a rectangular grid such that every point [x, y]  in the disk of radius rmax  in the x  direction and dy/2  in the y  direction of a grid point [x[i], y[j]] .  If [x[0], y[0]]  is the actual minimum (so in particular F(x[0],y[0]) <= CurrMin ), how much bigger than CurrMin  can F(x[i],y[j])  be?  Well, from Taylor's Theorem (remembering that the first derivatives of F(x,y)  at the minimum are 0), F(x[i],y[j]) = F(x[0],y[0])+F[x,x](x,y)*(x[i]-x[0])^2/2+F[x,y](x,y)*(x[i]-x[0])*(y[j]-y[0])+F[y,y](x,y)*(y[j]-y[0])^2/2
F(x[i],y[j]) = F(x[0],y[0])+F[x,x](x,y)*(x[i]-x[0])^2/2+F[x,y](x,y)*(x[i]-x[0])*(y[j]-y[0])+F[y,y](x,y)*(y[j]-y[0])^2/2  for some point x, y  on the line segment from [x[i], y[j]]  to [x[0], y[0]] .  In particular if F[x,x] <= a , abs(F[x,y]) <= b  and F[y,y] <= c , we have F(x[i],y[j]) <= CurrMin+a*dx^2/8+b*dx*dy/4+c*dy^2/8 .   

>    a:= op([1,2],evalf(evalr(subs(x=INTERVAL(-rmax..rmax),y=INTERVAL(-rmax..rmax),diff(F(x,y),x$2)))));

a := 18661.90914

>    b:= op([1,2],evalf(evalr(subs(x=INTERVAL(-rmax..rmax),y=INTERVAL(-rmax..rmax),abs(diff(F(x,y),x,y))))));

b := 100.

>    c:= op([1,2],evalf(evalr(subs(x=INTERVAL(-rmax..rmax),y=INTERVAL(-rmax..rmax),diff(F(x,y),y,y)))));

c := 144771.3802

With dx = .006 and dy = .002, I get

>    dz:=eval(a*dx^2/8 + b*dx*dy/4 + c*dy^2/8,{dx=.006,dy=.002});

dz := .1566642812

>    ti:= time():
dx:= .006: dy:= .002:
candidates:= NULL;
imax:= ceil(rmax/dx);
for i from -imax to imax do
  xi:= i*dx;
  if abs(xi) > rmax then jmax:= 1
  else jmax:= ceil(sqrt(rmax^2-xi^2)/dy)
  fi;
  for j from -jmax to jmax do
    yj:= j*dy;
    Fij:= evalhf(F(xi,yj));
    if Fij < CurrMin + dz then
      candidates:= candidates,[xi,yj,Fij];
      if Fij < CurrMin then
        CurrMin:= Fij;
        bestcand:= [xi,yj];
      fi
    fi
  od
od:
(time()-ti)*seconds;

>   

candidates := NULL

imax := 301

207.631*seconds

>    nops([candidates]);

75

>    CurrMin;

-3.30453521949094364

We've come down quite a bit from the previous CurrMin.  Now we can search closer to each of the 75 candidates.

>    dx:= .0006; dy:= .0002; dz:= a*dx^2/8 + b*dx*dy/4 + c*dy^2/8;candidates2:= NULL:

dx := .6e-3

dy := .2e-3

dz := .1566642812e-2

>    ti:= time():
for cand from 1 to 75 do
  for i from -5 to 4 do
    xi:= candidates[cand][1]+(i+1/2)*dx;
    for j from -5 to 4 do
      yj:= candidates[cand][2]+(j+1/2)*dy;
      Fij:= evalhf(F(xi,yj));
    if Fij < CurrMin + dz then
      candidates2:= candidates2,[xi,yj,Fij];
      if Fij < CurrMin then
        CurrMin:= Fij;
        bestcand:= [xi,yj];
      fi
    fi
  od
od
od:
(time()-ti)*seconds;
    

3.721*seconds

>    nops([candidates2]);

16

>    candidates2;

[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...
[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...
[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...
[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...
[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...
[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...
[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...
[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...
[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...
[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...
[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...
[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...
[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...
[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...
[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...
[-.2550000000e-1, .2105000000, -3.30321493672779853], [-.2550000000e-1, .2107000000, -3.30326071544024336], [-.2490000000e-1, .2099000000, -3.30356369591246635], [-.2490000000e-1, .2101000000, -3.30479...

>    plots[display]({
 plots[implicitplot](diff(F(x,y),x),x=-.0255-dx/2 .. -.0237+dx/2, y=0.2099-dy/2 .. 0.2111+dy/2, colour=red), plots[implicitplot](diff(F(x,y),y),x=-.0255-dx/2 .. -.0237+dx/2, y=0.2099-dy/2 .. 0.2111+dy/2, colour=blue) });

[Maple Plot]

>    Digits:= 30; fsolve({diff(F(x,y),x), diff(F(x,y),y)},{x,y},x=-.0255-dx/2 .. -.0237+dx/2,
y=0.2099-dy/2 .. 0.2111+dy/2);

Digits := 30

{y = .210612427155355770591591100555, x = -.244030796943751719036133083297e-1}
{y = .210612427155355770591591100555, x = -.244030796943751719036133083297e-1}

>    eval(F(x,y),%);

-3.30686864747523728007611377089

All 30 digits here are correct.