Help with Floating-Point Arithmetic

1 view (last 30 days)
Nick
Nick on 14 Oct 2011
I need to evaluate the expression,
X=(-0.5.*(1-A)+0.5.*((1-A).^2+(A*10).^2).^0.5).^0.5;
for an array of values, A, that ranges from about 1e-9 to 1e-8. Several of the results in X are 0 and I believe this is a result of A being swamped by 1. I'm trying to rearrange the arithmetic like in support note 1108, but I haven't had any luck yet.
Does anyone have a suggestion for getting past the round-off error?
Thanks.

Answers (2)

Walter Roberson
Walter Roberson on 14 Oct 2011
I would suggest tayloring it near 1e-9, but doing that in a symbolic engine with the digits set fairly high. Then convert it to variable precision arithmetic with a high number of digits, and rewrite that as a horner polynomial. Then drop the number of digits to get something you can evaluate at the MATLAB level.
If I bothered to do some numeric analysis, I might be able to figure out how many terms you needed to take it to, but midnight seems to have snuck up on me again.
  1 Comment
Nick
Nick on 18 Oct 2011
Thanks for your reply.
Let me see if I've got this straight, I would implement your steps as:
syms A
digits(2^10)
X=(-0.5*(1-A)+0.5*((1-A)^2+(A*10)^2)^0.5)^0.5;
T=taylor(X,1e-9);
V=vpa(T,2^10);
H=horner(V);
A=someValue;
digits(2^6)
eval(H)

Sign in to comment.


Jan
Jan on 14 Oct 2011
Reordering the terms to get a more stable calculation:
X = sqrt(-0.5 + 0.5 * A + 0.5 .* sqrt(1 - 2 * A + A .* A * 101))
The term (1 - 2 * A + A .* A * 101) remains critical if A is 1e-9, because the A.*A vanishes in the sum due to rounding.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!