How to solve this weird problem in Matlab?

Dear all,

I have a linear equation: a*x1-b*x2+x3=1.0. a and b are known parameters: a=b=10E+6. Now we know the solutions of x1 and x2: x1=1.001 and x2=0.9998. We want to solve x3.

The true solution of this equation is: x1=x2=x3=1.0. But because there are noise in the solution of x1 and x2, we got a weird solution of x3: x3=-1199.0, a very big value.

I am looking for a correct approach to solve for x3=1.

Thanks a lot. Bei

Answer 1

[John D'Errico]

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)
x3 =
       -11998.9999999981

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:

a*x1
ans =
   9980000    10020000
b*x1
ans =
   9980000    10020000

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)]))
x3 =
           -39999       40001

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:

x3 = 1 - (a*x1 - b*x2)

Then the mean of x3 will be 1, but the variance of x3 will be:

varx3 = a^2*(sigmax1^2 + sigmax2^2)
varx3 =
   200000000

sigmax3 = sqrt(varx3)
sigmax3 =
         14142.135623731

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.

Answer 2

[YT]

So I've tried this and got the same result.

syms x3
a = 10e6;
b = 10e6;
x1 = 1.001; 
x2 = 0.9998;
eqn = a*x1 - b*x2 + x3 == 1;
x3 = vpa(solve(eqn,x3));
%-11998.999999998137354850769042969

I was confused at first, but if you substitute this value into the equation for x3, it will result in 1=1, so its `solve` did a good job.

After this I started to examine the function a bit better and the values you're trying to put in (deviation in x1 and x2; large values for a and b). So instead of a=b=10e6, let's take a=b=1, then the result for x3 would be

0.99880000000000013216094885137863

Now I'm not an expert at making conclusions when it comes to math, but I would say that if your looking for the outcome x3 = 1 in the equation a*x1 - b*x2 + x3 == 1, you have to make sure that x1 = x2. Even the slightest deviation between x1 and x2 won't get you x3 = 1.

[Benson Gou]

We actually do not know the exact solution of x1 and x2. This example was created intentionally to describe the problem.

Thanks. Benson