Linsolve broken recently? was Re: Oscillatory integrals

On 09/27/2013 07:31 PM, Richard Fateman wrote:
>  This new version makes LevinInt work:

Works for me now. I have been listing a few integrals to
test with different methods before coding. One example is


integrate(cos(x^2),x,0,inf) or integrate(cos(x^2),x,0,1000)

I tried changing variables and rescaling the interval etc. to get
different routines to work.  I could not get LevinInt to perform well
on this problem.  quad_qagi can't include the origin. quad_qago is a
bit better than quad_qags, but also can't handle a singularity. But
quad_qags can include the origin and integrate over a very large
interval.

But I haven't compared the more complicated things like cos(sinh(x))....
Maybe there are some that quadpack cannot do well, but LevinInt or
something like it can.

Then there are the other routines.... I think this is a big
project.

I wonder if changevar has always been broken in that it gives
the wrong sign sometimes. Or is that sort of recent? There
are a couple of bug reports on it.

>
>  LevinInt(f,g,a,b,omega,k):= 
block([varlist,p,pointlist,sol,ratprint:false],
>  varlist: makelist(alpha[i],i,0,k-1),
>  p: varlist . makelist(x^i,i,0,k-1),
>  pointlist : makelist(a+ (b-a)*i/(k-1),i,0,k-1),
>  define(eq(x), diff(p,x)+%i*omega*diff(f,x)*p-g),
>  sol:linsolve(eqs:map(eq,pointlist),varlist),
>  define(H(x), exp(%i*omega*f)*float(expand(ev(p,sol)))),
>  rectform(H(b)-H(a)))$
>
>
>  The only change is from
>
>  define(H(x), exp(%i*omega*f)* ev(p,sol)),
>
>  to
>  define(H(x), exp(%i*omega*f)*float(expand(ev(p,sol)))),
>
>  strictly speaking, the float() or it could be bfloat() is not necessary,
>  but the expand() is. I think linsolve provides bad (ill-conditioned) 
pivot elements when
>  coefficients are complex.
>  ...
>
> > (%i3) LevinInt(sinh(x),cos(x),-1.0,1.0,1000,7);
>
>  The results from linsolve look substantially different from Maxima 
5.25 and 5.29.
>  They apparently handle this somewhat ill-conditioned system in different
>  ways. I think the new version is not as good as the old one on this
>  system, but I'm not sure why. Shouldn't the exact solution to an exact 
system
>  be the same?
>  Anyway, the new version seems to give rectform indigestion.
>
>  The new version may produce a kind of factored form that is like
>  tiny*(huge+ huge*%i)
>
>  I'm hoping that someone else will diagnose this better.
>
>
>