Re: Bracketing Issues with fzero

[Jaroslav Hajek]

On Fri, Mar 26, 2010 at 10:19 PM, Fabio Somenzi <address@hidden> wrote:
> I tried to work around the fzero/fsolve problem by supplying the
> Jacobian to fsolve.  I found that my installation of 3.2.4 did not
> even run the example in the documentation of fsolve. It complained
> about incorrect "cells." I uninstalled 3.2.4, installed 3.0.5, and all
> is well.
>
> This may be some weirdness specific to the Windows port, but I thought
> I'd report my findings.  Speaking of which, in my first mail I forgot
> to mention that with 3.2.4, when fsolve returns -3, perror complains
> of an unexpected value of info.
>
> Fabio
>
> On Fri, Mar 26, 2010 at 10:33 AM, Fabio Somenzi <address@hidden> wrote:
>> Hi,
>>
>> I just installed the stand-alone version of Octave 3.2.4 on my netbook
>> running Windows 7 Enterprise.
>>
>> The script below, with fsolve instead of fzero, runs fine on Octave
>> 3.0.3 which I have installed on other machines.  With 3.2.4, fsolve
>> returns -3 in many calls.  After consulting the documentation, I
>> switched to fzero, but I get messages like
>>
>> fzero: zero point is not bracketed
>> fzero: not a valid initial bracketing
>>
>> I thought I had a work-around when I realized that the first message
>> (zero point is not bracketed) occurred when the guess for the solution
>> was in fact the solution.  (This happens every time at the first
>> iteration of the inner loop).  So, I thought I could just slightly
>> perturb the initial value or skip the first (needless) iteration of
>> the inner loop, but I only end up getting the second message (not a
>> valid initial bracketing) later on.
>>
>> How should I proceed?  This is my first post.  Feel free to advise on
>> how to report problems effectively.
>>
>> Thanks
>>
>> Fabio
>>
>> -------------------------------------------------------------------------------------------------------------
>> ## Model for a falling climber in which we assume friction at the top
>> ## biner.  This is a simplified model based on conservation of energy.
>> ## This model uses the length of unstretched rope that slips through the
>> ## carabiner as basic variable.
>>
>> ## The fall distance is assumed to be twice the length of rope initially
>> ## on the climber's side of the carabiner.
>>
>> global m = 80;          # kg    : mass (80)
>> global g = 9.81;        # m/s^2 : acceleration due to gravity
>> global k = 20000;       # N     : rope modulus
>> global mu;              #       : 1 - pulley efficiency of the carabiner
>> global L1;              # m     : length of rope on the climber's side
>> global L2;              # m     : length of rope on the belayer's side
>>
>> ## Estimate impact force according to Wexler's equation.
>> function iforce = impactf(ff)
>>
>>  global m g k;
>>
>>  w = m * g;
>>  iforce = w * (1 + sqrt(1 + 2 * k * ff / w));
>> endfunction
>>
>> ## Energy balance equation.  For given slippage s, it returns the energy
>> ## imbalance.  We seek the value of s that makes that imbalance zero.
>> function w = f(s)
>>
>>  global m g k mu;
>>  global L1 L2;
>>
>>  z = L2 - s;
>>  q = 1 - mu;
>>
>>  w = k/(2*q^2) * (mu*(2-mu)*s + (L1 + L2*q))*(s/z)^2 - \
>>      m*g*(2*L1 + (L1+s)*s/(q*z) + s);
>>
>> endfunction
>>
>> ## Setup.
>> muvals = [0, 0.25, 0.5, 0.75, 0.999]; # values of mu to be considered
>> L1 = 2; # To vary the fall factor, we keep L1 fixed and change L2
>> n = 400; # Number of samples of fall factor
>> ff = linspace(0.00001,1.99999,n);
>>
>> ## Computation of one curve for each value of mu.
>> for i = 1:n
>>  L2 = (2/ff(i) - 1) * L1;
>>  rope(i) = L1+L2; # length of rope required by the given fall factor
>>  ## To avoid numerical problems, we need a reasonable initial value. We
>>  ## use the analytical solution to derive one for mu = 0 and then we
>>  ## use the solution for one value of mu as initial value of s for the
>>  ## next value of mu.  This works best if the values of mu are
>>  ## considered in increasing order.
>>  nof(i) = impactf(ff(i));  # Rope tension according to Wexler's equation
>>  guess = L2 * nof(i) / (nof(i)+k); # corresponding slip
>>  for j = 1:length(muvals)
>>    mu = muvals(j);
>>    [s,value,info,myoutput] = fzero(@f, guess);
>>    if (info != 1)
>>      if (info < -4 || info > 1)
>>        printf("fzero returned %d\n", info);
>>      else
>>        perror("fzero", info);
>>      endif
>>    endif
>>    guess = s;
>>    slip(i,j) = s/L2; # slippage
>>    T2(i,j) = k * s / (L2-s); # force on belay
>>    T1(i,j) = T2(i,j) / (1-mu); # force on climber
>>    T(i,j) = T1(i,j) + T2(i,j); # force at the anchor
>>    stretch(i,j) = ((L1+s)*s/((1-mu)*(L2-s)) + s)/(L1+L2); # total stretch
>>  endfor
>> endfor
>>
>> ## Plotting of curves.
>> plot(ff, T, "linewidth", 2);
>> for j = 1:length(muvals)
>>  strmat{j} = strcat("mu = ", num2str(muvals(j)));
>> endfor
>> legend(strmat, "location", "northwest");
>> #print('topAnchor.eps', '-depsc');
>>
>


I see fzero was breaking down in the (rare) case when the bracketing
included an exact root; got that fixed. I also improved the docstring
for fzero so that it explains that a bracketing is expected (i.e. both
a below-zero and above-zero point).
If you supply just a single point (which you do in your script), fzero
tries to determine the bracketing, but if it fails, it ends in an
error because the algorithm can't continue.
For most functions, bracketing the root is trivial. If it can't be
done (function is almost always positive or always negative), consider
using fminbnd instead and search for a minimum/maximum.

I should also note that your f function appears to be a cubic
polynomial in s; hence, it may be a much better idea to find the root
using "roots".

hth


-- 
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz