How can I find local and global minima of a 2D non continuous function

Hello,

I am currently trying to compute both local and global minimum of a two dimensional function. I already tried several methods but it seems that the function's definition I am using cannot be used here. The different MATLAB functions I use like imregionalmin, fminunc, GlobalSearch etc... return me wrong results (multiple global minima) that vary depending on the mesh size I use, or the starting point of the function. You can find below the function itself. I know from references and experiment that the result should be a unique global minima and one local minima. if true

    b = pi;
    A11 = 6.91;
    D11 = 0.081;
    D12 = 0.037;
    D22 = 0.081;
    R = 15;

    fonc = @(x,y) (0.5*b*R).*(D11.*x.^2 + 2*D12.*x.*(y - 1/R) + D22.*(y - 1/R).^2)...
    + (0.5*A11).*(0.5*b*R.*(x./y).^2 + 0.5*sin(b*R*y)*(x.^2./y.^3) -...
    4*(sin(0.5*b*R*y)).^2*(1/(b*R))*(x.^2./y.^4));

end

[Walter Roberson]

My tests indicate that there are minima at

[0, 1/15] -- global minima

[0.65583738972311451561151771866332878854679817998043e-2, 0.2438773518742782030725331e-2]

[0.30235293603562578703919586060332531258654382992792e-1, 0.98750479772327783826276243144452399673447717922487e-4]

[0.30449374732803001519898620416944969463601022607734e-1, 5.8786799832110556917894902694773842319226264776278*10^(-6)]

[0.30452674896762953857019895149460257634111447217088e-1, 7.8020009110210131711531349506573007830550822287214*10^(-13)]

To find this, I used symbolic x and y, and took

bestx1 = solve(diff(fonc,x),x);
F = subs(fonc(x,y), x, bestx1);
eqn = diff(F,y);

Then if you look at the zeros of eqn very closely, you get the y values from the above sets, and you can substitute those into bestx1 to get the corresponding x.

Now... an important point along the way is that this is for the fonc as you have defined you. You have defined it using floating point numbers such as pi (= 884279719003555/281474976710656 in internal representation) and 6.91 (= 1944992089070633/281474976710656 in internal representation) . pi in your formula is not exactly the same as the transcendental constant π and 6.91 is not exactly the same as 691/100 .

If you redefine the floating point constants as-if they were integers divided by powers of 10, but making pi the transcendental constant, then you get a different equation that is tougher to find formulas for to narrow down the zeros.

[Walter Roberson]

I ran vpasolve() with a number of starting points. It confirms values of approximately 9.87705604732379e-05 0.00243877351724839 for y but does not confirm the 5E-6 and 7E-13 values that I found with other software.

Solving only appears to work nicely in one direction.

    Q = @(v) sym(v, 'r');
    b = sym('pi');
    A11 = Q(6.91);
    D11 = Q(0.081);
    D12 = Q(0.037);
    D22 = Q(0.081);
    R = 15;
    H = Q(0.5);
    fonc = @(x,y) (H*b*R).*(D11.*x.^2 + 2*D12.*x.*(y - 1/R) + D22.*(y - 1/R).^2)...
    + (H*A11).*(H*b*R.*(x./y).^2 + H*sin(b*R*y)*(x.^2./y.^3) -...
    4*(sin(H*b*R*y)).^2*(1/(b*R))*(x.^2./y.^4));
    syms x y
    F = fonc(x,y);
    bestx1 = solve(diff(F,x), x);
    FF1 = subs(F, x, bestx1);
    besty1 = vpasolve(diff(FF1,y), y)
    besty2 = solve(diff(F,y), y);
    FF2 = subs(F, y, besty2);
    bestx2 = vpasolve(diff(FF2,x), x)

    Y = sym(logspace(-20,3,100));
    for K = 1 : length(Y); sol(K) = vpasolve(dF1, y, Y(K)); end
    uniquetol( double( [sol{:}] ) )

Answer 1

[Zakariya LOULOU]

Thanks for your answer !

I understand what you mean by how I should define my numbers. I guess the best thing to do right now is to stick with floating point numbers to solve fonc. Using your method, I found the global minima (0, 1/15), which is relevant with the experiment I did. However, in order to fund the zeros of eqn, I am trying to use fzero function, with a lot of differents starting points located within (-0.05; 0.25) (interval of interest here) and I don't find the same results as yours.