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

dlakelan wrote:
> 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
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
> .
>
>   
Don't know about your functions since you have not given them explictly, 
but I have noticed that putting complicated functions into one big 
procedure has caused problems for me.

What seems to work is to use temporary variables to take the result of 
intermediate computations, and then assign them to later functions. 

To be overly simplistic,  suppose that we have the composition  h(x):= 
f1(f2(x)) of two functions (in reality the routines may be much more 
complicated of course).

If f1 takes a short time to execute and f2 takes a long time to execute 
(or either depends on free or bound variables in some way only maxima 
deciphers), then there may be errors produced.  This happens a lot when 
one of f1 or f2 has some 'makelist' routines embedded inside. So, my 
solution would be to explicitly separate the routines (or functions ) as:

 g(x):= block([temp_1,temp_2], temp_1: f2(x),  temp_2: f1(temp_1), 
return(temp_2))$

I have encountered situations where  the simple routine
h(x):= f1(f2(x)) 

would cause errors, but the equivalent routine g(x) above would work.

Another comment:
Plot2d is especially sensitive to many things.  Why not try draw2d instead?

HTH
-sen