If your aim is to get x == X, where x denotes numerical solution of AA*x=B for AA and B defined in your example, then it is impossible.

Firstly some notation, which I'll use:

Ax=b; A=[aij]; x=[xj]; b=[bi];

s( ) = absolute error

r( ) = relative error

A\b = arbitrary method solvin linear system

~ means "similar" or "comparable to", e.g. x~y is meant as x is of the similar order as y

AA, B, X = variables defined in your example

Obtained x=AA\B (x~=X) is correct solution at the used computational precision. It can be shown by

x = AA\B; sB = AA*x - B,

where sB is typically in order of 1e18 (note that result of (1e34+5e17)==1e34 is true).

In case of double the precision should be 16 digits, so for the example round(log10(A(i,:)) == [34, 27, 19, 12, 4] and round(log10(bi)) = 34, Ax = b+s(b) gives s(b)~1e18 and any aij*xj<s(b) becomes insignificant.

It restricts also the number of valid digits for all aij and xj. If error level (variance of Ax and/or b) of Ax=b is given as s(b)^2 = sum_j{[aij*s(xj)]^2+[s(aij)*xj]^2} and numbers are represented by a decadic decomposition as y=sum_k(y(k)*10^n(k)), its obvious that any [aij(k)*10^n(k)].*[xj(q)*10^n(q)] < 1e18 turns to be negligible in sense of equality Ax==b.

For orders of magnitude of aij [1e34,1e27,1e19,1e12,1e4] the highest obtainable precision on x is [1e-16,1e-9,1e-1,1e6,1e12]. It's clear that x(4) and x(5) could be never estimated correctly, if expected value is in between 0 and 1.

For last precise digit of ai1 with order of 1e34/1e16=1e18 the highest precision applicable to x1 is 1e0, thus even for the coefficient of the highest order not all valid digits of the coefficient are significant. For requested precision of x1 say 1e-4 the last significant digit of ai1 is at place of about 1e22.

The above conclusions about obtainable precision of x can be illustrated by a single step of Jacobi's point-relaxation iterative method with X as the initial guess.

D = diag(diag(AA));

xe = inv(D)*((D-AA)*X + B);

The difference between X and xe is in perfect line with the maximal expectable precisions of [1e-16,1e-9,1e-1,1e6,1e12].

The first idea how to reach the aim x==X is to use higher precision of calculation. It can be found some tools to do that (e.g. <a href="

http://gmplib.org/">

http://gmplib.org/), but since I'm not familiar with those and it was particularly what attracted my attantion I wrote something myself (see uploaded files below). Calculation with variable precision will be noted infP(). Surprisingly x=infP(AA\B) is close to AA\B instead of X (note infP(A\b) uses simple Gauss's elimination). The reason is in the fact that B generated by AA*X is already loaded with the error of order 1e18 as can be seen by comparison of B=AA*X and infP(B=AA*X). Solving the system x=infP(AA\infP(B=AA*X)) gives the desired result x==X (in double prec.).

The conclusion can be made that the solution x=AA\B is correct for given computational precision, so that A*x = b + s(b) with the maximal precision (in sence of orders) of s(xj) = log10(aij)/s(bi). To be sure that all estimated xj are grater than their s(xj), we can quess the needed calculation precision by demanding the lowest acceptable precision of x to be grater or equal to the round-off error of calculation, i.e. s(xn)>=max_j(s(aij*xj))~s(bi), where n denotes the index of the component with coefficient of the lowest order (in example n=5 and s(x5)>=s(ai1*x1)~s(bi)). It can be expressed as

(x1*s(a1))^2 + (s(x1)*a1)^2 <= (xn*s(an))^2 + (s(xn)*an)^2

(x1*ar*r())^2 + (x1*a1*r())^2 <= (xn*an*r())^2 + (s(xn)*an)^2

r() = r(aj) = r(xj) = r(b) = relative precision of calculation, e.g. for double r()=1e-16

it could be assumed that x1*ar*r() == xn*an*r(), then

s(xn) is left unchanged, since it represents the lowest requested precision and not a round-off error

given by computation

x1*a1*r() <= s(xn)*an

x1*r() = s(x1) <= s(xn) * an/a1

if we assume x1~1, then r() <= s(xn) * an/a1

So if we want x5 to be precise up to 1e-4, we heve to use r() <= 1e-4 * 1e4/1e34 = 1e-34. It is only rough guess, therefore smaller r() is desirable. The important point is that with this precision it must be estimated/calculated/measured/known also the righthand side (b) as shown above. I feel this condition as the major obstacle in practical problems. If it is not your case, then I hope it helps you.

File with code for variable prec. calculations:

operation.m It isn't "all input errors proof" and also it isn't optimized anyhow, so please be patient during evaluation.

I wrote it because of curiosity and not for a real use. It's advisable to use some existing package or tool for variable precision calculations in case of solving your real problems.

File with demo / checking of operation.m performance / correctnes:

OperPerformance.m Note the last example of estimation of e, how the conversion to double changes a precise number.

File dealing with your example:

AAx_B_demo.m The order of printouts follows text of this reply. Actualy this is the only script you need to run (assuming the operation.m is saved on the path).