Fastest way to minimize a function

I have to perform many minimizations of a function, so I want to do it quickly. Here is the problem: given different values of T, the minimization should find L = ((g*T^2)/(2*pi))*tanh(2*pi*D/L) where g and D are constants. Below is how I am doing it, but is there a faster way?

g = 9.8066 ;
D = 50 ; 
periodList = linspace(3, 15, 100) ;
fh = @(L,T) (L-((g*T.^2)/(2*pi))*tanh(2*pi*D/L)) ; % Handle to L as a function of T
for iPeriod = 1:length(periodList)
    T = periodList(iPeriod) ; 
    guessL = g*(T^2)/(2*pi) ; % Initial guess at L
    L(iPeriod) = fzero(@(L) fh(L, T), guessL) ;
end

Since this problem is finding the wavelength L of water waves given the wave period T and water depth D, there are some constraints on the problem: L must always be positive with a nominal solution accuracy of +/-1 meter, the period T will lie between 0 and ~20 seconds, and the water depth D ranges from ~10 to 1000 meters. g is the acceleration due to gravity, so it will always be 9.8066 m/s^2.

[Sean de Wolski]

fmincon() or ga() could both be your friends.

[K E]

Are these faster than fzero (just curious since I don't have the Optimization Toolbox)?

[Sean de Wolski]

ga is definitely not faster, fmincon might be but probably not. But they're significantly more powerful.

Answer 1

[Walter Roberson]

Unfortunately the solutions to your function suffer from sufficient arithmetic problems as to make minimization not of value: there are upper and lower bounds beyond which all L values produce the same result (different for the upper and lower bound though.)

The number of real roots of L vary considerably; with the test values I was using for g and D, I found up to 45 real roots. Unfortunately due to precision issues, the roots were indistinguishable except that one of them was the positive value corresponding to the other (negative) values.

With the constants you show above, all of the really interesting behavior occurs between T=0 and T=2, with some semi-interesting behavior up to T=5, and boring behavior after that. The positive L values increase as T increases -- and contrawise, the negative L values corresponding get more negative, so if fsolve() happens to hit the negative root, you will find the numbers becoming more negative (more minimum) as T increases. If you want the T that has the minimum L, then you might as well just choose the upper bound of your T range once T has gotten large enough that your results are dominated by floating point round-off. By the way, if you work matters through from the point of view of round-off then the best initial guess is the exact negative of the guess calculation you show.

I am having my system chase some obscure behavior near T=2; I do not know yet whether it will be meaningful or not.

[Walter Roberson]

The behavior value near T=2 was a red herring, a glitch in the system I was using.

[K E]

I will work through the useful info in your answer, but meanwhile there are constraints on the solution (such as L>0) that I did not mention originally. I have added those constraints to the original question. Perhaps they eliminate some of the regions in which the solution behavior is odd; I'll check.

[Walter Roberson]

The form of the solution turns out to simplify a fair bit and be more arithmetically stable if you substitute L -> g*(T^2)/(2*pi) + dL and solve for dL -- solve for the difference from the guess. Most of the roots for dL will still be strongly negative, so you have to pick out the values for which g*(T^2)/(2*pi) + dL >= 0, and then pick the smallest root out of that list.

Using this method and the sample parameters of your Question, dL heads towards -200-ish near t=20, and does so faster than exponential. Unfortunately the version I have now still is not quick, probably because I am having Maple find the entire set of roots and picking through it, when it would probably be faster to ask for a numeric solution near dL = 0.

[Walter Roberson]

Bleh, the number of roots of the equation depends upon the number of digits of precision you are using.

It appears at the moment (but no proof!) that when an appropriate solution exists, it is based upon the only positive solution for the RootOf(), and will be difficult to find numerically unless you constrain the boundaries of solution to be positive (the function is very steep at the positive root, so solving by way of a guess will usually miss the positive value.) On the other hand, it will also be the negative of the easily-found negative root of the RootOf()

[K E]

Thanks so much for your insights. Again there was information that I should have included that avoids the problem. When the wave period gets very long (t>20s) and L becomes negative which is unphysical, then you can assume the wave obeys the approximation L=sqrt(g*D)*T which applies when the water is shallow relative to the wave's dimensions.

[Walter Roberson]

The missed roots show up most clearly in about the area of T=Pi .

The number of roots identified varies wildly with very small changes in T, from 1 up to (at least) 72. Nasty precision problems.