>>>>> "Robert" == Robert Dodier <robert.dodier at gmail.com> writes:
Robert> On 2013-08-28, Sean Lake <odysseus9672 at gmail.com> wrote:
>> f(x,y) := exp( - ( y - x - log( 1 +exp(y))/3)^2 / 2);
>> N(x) := quad_qagi( f(x, y), y, minf, inf )[1];
>> plot2d( N(x), [x, -5, 5] );
>>
>> When I try to evaluate N(0.0), I get what looks like an unevaluated
>> noun form
Robert> quad_qagi tries to evaluate log(exp(y) + 1) for large y, and gets a
Robert> floating point error, at which point it gives up. (Try :lisp (trace
Robert> simpexpt) and then call quad_qagi to see that.) The only way around
I think it would be nice if maxima indicated that somehow instead of
just returning a noun form. Users should have to trace internal
functions to figure that out.
Robert> it is to rephrase the integrand so it doesn't cause an error. Maybe
Robert> you can replace log(exp(y) + 1) by a rational approximation or
Robert> just y for large enough y.
Robert> Incidentally the part from minf to 0 doesn't seem to cause any trouble.
A plot of f(0,y) shows that the f(0,y) is pretty much concentrated
between -5 and 5. And f(x,y) is roughly f(0,y) translated by x, so
the infinite integral can be replaced with a finite one.
quad_qags(f(0,y),y,-10,10) is probably a very good approximation to
the quad_qagi(f(0,y),y,minf,inf). You'll have to do the analysis to
figure out how good the approximation is.
Ray