[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Bracketing Issues with fzero
From: |
Fabio Somenzi |
Subject: |
Re: Bracketing Issues with fzero |
Date: |
Fri, 26 Mar 2010 15:19:53 -0600 |
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');
>