Argh, why can't I ever figure out plot2d/find_root/etc



Stavros Macrakis wrote:
> On Thu, Mar 6, 2008 at 11:05 PM, dlakelan <dlakelan at street-artists.org> wrote:

>>  Why does this suck so much?
> 
> Despite your giving a useless bug report with no details, and using
> vulgar language to insult everyone who's worked on Maxima over the
> years, I suspect you will actually get help on this list.

Stavros:

I am sorry to have offended you or anyone else. As you can imagine 
computers can be frustrating things at times. I am a very enthusiastic 
Maxima user, and I want things that look right to work because I want to 
be able to recommend maxima to others... But I run into this kind of 
confusion almost every time I try to use maxima for this sort of thing.

The simple example is not very easy to give, because these things 
USUALLY WORK for the simple case... Kudos to the maxima developers for that!

Here's a description of what i'm trying to do. perhaps someone will be 
able to come up with a simple case that shows the bug.

I have a model for concrete stress/strain relationship. I am analyzing 
the bending properties of a concrete beam. I apply a degree of curvature 
to the beam, and then I try to balance the tension force in the steel 
with the integral of the compressive stress distribution over the cross 
sectional surface of the beam. The strain is assumed to be linear with 
distance from the "neutral axis". So I need to find the neutral axis...

neutral_axis(phi,w,d) := find_root('(totalforce(a,phi,w,d)),a,0.0,d);

w, and d are the width and depth of the rectangular cross section... phi 
is the angle of curvature

As you can see I've already started quoting things to try to debug 
what's going on.... but it hasn't changed the behavior...

totalforce is a function which calls quad_qag to integrate the stress 
over the rectangular cross section given a neutral axis depth at "a"...

Moment is a function which calculates the total internal moment, which 
is the torque around the neutral axis applied to one side of the beam by 
the other side. It calls neutral axis to find the neutral axis, and then 
it calls quad_qag to integrate the torque caused by the stress 
distribution in the concrete, and adds the torque caused by the 
reinforcement steel.

When I call Moment at the command prompt, it produces a nice floating 
point number. when I call

plot2d(Moment(phi,w,d),[phi,.0001,.002]);

I get an error about an expression involving "$A" not being a floating 
point number, which makes me think that find_root deep inside 
neutral_axis is unhappy...

But when I create a list of x values and makelist(Moment(i,w,d),i,xvals) 
I can plot the discrete x,y points no problem.

When I start debugging this sometimes I can change what the error is... 
sometimes it will be an error about the stress distribution function 
(which is a function of the form
if(this) then that else if (the other) then formula;

This suggests to me that when this alternative error occurs, quad_qag is 
complaining about its first argument not evaluating to a float, similar 
to find_root.

In all cases where I get an error when plotting, or when doing numerical 
integration, etc I get a perfectly well formed floating point number 
when i call the function at the command line...

I can imagine that we can emulate this behavior by simplifying down, so 
i will try to post a simplified version that still shows the behavior.

In the mean time, I originally thought that this problem might be common 
enough that someone would say "oh yes, you need to be careful that your 
argument is quoted in a certain way, or use funmake or..." which is why 
I stupidly posted my rash comment.

Thanks to the whole maxima team for making a great piece of software, 
here's hoping your QA department (me) will be more helpful in the future.

Dan