quad_qagp behavior
- Subject: quad_qagp behavior
- From: Edwin Woollett
- Date: Sat, 21 Jan 2012 14:12:45 -0800
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))
(%i8) for y in [1,2,3] do
print(y,float(r1(y)),float(r2(y)))$
1 1.090135361218102 1.090135361218102
2 1.048689739767972 1.048689739767972
3 2.661985611481998 2.661985611481998
-------------------------------------------
Ted Woollett