Re: Bracketing Issues with fzero

[Fabio Somenzi]

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