On 1/21/12 2:12 PM, Edwin Woollett wrote:
>
> testing the new quadpack utility, quad_qagp, weird result not understood:
>
> e = realpart(1/e1) gives correct result
>
> e = 1/realpart(e1) gives really bad answer
> -----------------------------------------------------------------------
> Maxima 5.26.0 http://maxima.sourceforge.net
> using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (a.k.a. GCL)
>
> ok:
> (%i1) quad_qagp (realpart (1/sqrt (sin(x))), x, 1, 5, [%pi]);
>
> (%o1) [3.209309789377543,6.1241798299249695E-9,294,0]
>
> not ok:
> (%i2) quad_qagp (1/realpart (sqrt (sin(x))), x, 1, 5, [%pi]);
>
> (%o2) [4.7552920825719432E+16,7.0149712E+7,294,0]
>
>
> (%i3) realpart(1/sqrt(sin(x)));
>
> (%o3) cos(atan2(0,sin(x))/2)/sqrt(abs(sin(x)))
>
> (%i4) 1/realpart(sqrt(sin(x)));
>
> (%o4) 1/(cos(atan2(0,sin(x))/2)*sqrt(abs(sin(x))))
>
> compare the two resulting expressions numerically:
>
> (%i5) r1(x):= realpart(1/sqrt(sin(x)))$
>
> (%i6) r2(x) := 1/realpart(sqrt(sin(x)))$
>
> (%i7) for y in [1,2,3] do
> print(y,r1(y),r2(y))$
>
> 1 1/sqrt(sin(1)) 1/sqrt(sin(1)) 2 1/sqrt(sin(2)) 1/sqrt(sin(2)) 3
> 1/sqrt(sin(3)) 1/sqrt(sin(3))
You didn't take enough samples. I plotted r1(x) - r2(x) for x from %pi
to 5. and the plot is not zero. It is zero for x from 1 to %pi.
So I guess realpart is buggy in this case, although ratsimp(r1(x)-r2(x))
is zero.
Ray