A technical term for this is massive subtractive cancellation. I suppose you could also call it an ill-posed problem. So what is happening? Is there a way to resolve it?
x1 = 1.001;
x2 = 0.9998;
a = 10e6;
b = 10e6;
x3 = 1 - (a*x1 - b*x2)
Why would you be surprised? That you know that you want a solution is not really relevant. a and b are sufficiently large that the product of a*x1 and b*x2 are relatively big numbers. The uncertainty in x1 and x2 is sufficiently large that a and b magnify the noise.
Perhaps what we need to do is interval arithmetic here to understand the problem. Sadly, there are not really any interval arithmetic tools in MATLAB, but it is not difficult to implement a simple set of tools.
Thus, suppose we write x1 and x2 as the intervals
x1 = [0.998 1.002];
x2 = [0.998 1.002];
where x1 and x2 could be ANY number in the corresponding interval. I just picked a range for each of x1 and x2. Then we have a*x1 as:
Computing the difference of two intervals, what can we infer about x3?
x3 = sort(1 - ([a*x1(1) - b*x2(2),a*x1(2) - b*x2(1)]))
So x3 is SOME number in that region. That is as much as you can say.
You CANNOT infer that x3 is a number around 1.
Another way of looking at this is to assume that x1 and x2 are random variables, with a mean of 1, but lets say they have a standard deviation of roughly 0.001. So we don't really know what they are. We only have a mean and variance.
mux1 = 1;
mux2 = 1;
sigmax1 = 0.001;
sigmax2 = 0.001;
What is the mean and variance of the number x3, if we assume that x1 and x2 are normally distributed random variables with the indicated distribution parameters? Thus, if we have:
Then the mean of x3 will be 1, but the variance of x3 will be:
varx3 = a^2*(sigmax1^2 + sigmax2^2)
sigmax3 = sqrt(varx3)
Plus or minus 2 or 3 sigma there is pretty big.
The point is, the uncertainty in x1 and x2 is simply too large to resolve x3 with any degree of certainty.