Am Dienstag, den 09.12.2008, 23:15 -0700 schrieb Robert Dodier:
> On 12/9/08, Edwin Woollett <woollett at charter.net> wrote:
>
> > assume (t > 0);
> > x (t) := ''(integrate (cos (s^2 * %pi/2), s, 0, t));
> > y (t) := ''(integrate (sin (s^2 * %pi/2), s, 0, t));
>
> > Thank you for reminding me of this construct, which also works
> > with draw2d. However things fall apart if you go past about
> > t = 4.8, due to a current bug in the evaluation of erf(z).
>
> > (%i9) xfr(t) := realpart(float(x(t)))$
> > Evaluation took 0.0000 seconds (0.0000 elapsed)
> >
> > (%i10) map('xfr,[1,2,3,4,5,6,7]);
> >
> > (%o10)
> > [0.77989340037682,0.48825340607534,0.60572078930926,0.49842588395517,
> > 0.92294250203326,5394606.134146459,1.1458348613507996E+16]
>
> Hmm, this is very curious. It's not clear exactly what's happening in
> the erf function -- maybe x(t) is the difference of erf values which have
> very large magnitude, in which case inaccurate floating point results are
> not too surprising, or maybe erf is just incorrect for some arguments.
Yes, the erf function has a problem with some regions of complex values.
I am working on the problem as reported on this mailing list.
The error is in the gamma_incomplete function. The series expansion
converges to a wrong result in some parts of the complex plane. For
small arguments all is correct, but the erf(z) function need to
calculate the Incomplete Gamma function for an argument z^2. Therefore
the error in the series expansion increases very fast for the erf
function.
I am continuing the work to correct the numerical routine of
gamma_incomplete. Here a first result with the corrected code for the
above example:
(%i8) x(t):=realpart(float(''integrate(cos(s^2*%pi/2),s,0,t)))$
(%i9) map('x,[1,2,3,4,5,6,7]);
(%o9) [.7798934003768226, .4882534060753405, .6057207892976856,
.4984260330381772, .5636311887040121,
0.499531467855501, .5454670925469698]
Dieter Kaiser