How to avoid rounding error

When I use lu[A] for A = [10^(-20) 1 ; 1 2] I get 2 matrices (L and U). When I multiply them, the result is not the same as the original matrix A. What method can I use in order to get the correct matrix A?

[Rik]

You could try vpa.

The more fundamental problem is that computers have finite precision. If you want infinite precision, you will need to use algebraic tools. Not every problem can be solved perfectly. The general solution for this is to avoid problems that span more than 20 orders of magnitude, so you can rely on eps to estimate if your results are close enough.

Answer 1

[Walter Roberson]

A = [sym(10)^(-20) 1 ; 1 2]

A = 

[L,U] = lu(A)

L = 

U = 

L*U - A

ans = 

You can see from this that in order to get back A exactly, then you need a system that can distinguish 99999999999999999998 from 100000000000000000000, but

eps(100000000000000000000)

ans = 16384

it certainly is not double precision arithmetic.

[Walter Roberson]

syms N real
A = [sym(10)^(-N) 1 ; 1 2]

A = 

[L,U] = lu(A)

L = 

U = 

eqn = U(2,2) == -1/eps

eqn = 

solve(eqn)

ans = 

vpa(ans)

ans = 

15.653559774527022343979915836331

So beyond about 10^15.65 you go beyond what can be represented exactly in double precision.