Argh, why can't I ever figure out plot2d/find_root/etc
Subject: Argh, why can't I ever figure out plot2d/find_root/etc
From: dlakelan
Date: Fri, 07 Mar 2008 09:33:21 -0800
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