Bracketing Issues with fzero

[Fabio Somenzi]

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');