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.
I have looked into whether it's possible to rephrase x(t) and y(t) to
avoid computing the difference of erf's.
The best I can do is x(t) = (1/2) (Re foo - Im foo),
where foo = erf((1 - i) t sqrt(pi)/2), probably a similar expression for y(t).
But I can't figure out a relation between Re erf and Im erf
or an efficient way to evaluate the difference. Ideas?
best
Robert Dodier