Effect of Numerical Precision on Solving Linear Equations

Also, introduction to Mathematica!

Setting up the equations.

We define a system of linear equations in Mathematica. We are interested in the system:

0.1036 x + 0.2122 y = 0.7381
0.2081 x + 0.4247 y = 0.9327

which is given by the intersection of two lines in the plane. This is a matrix equation in the form 

Ax = b

Mathematica will solve this (and many much larger equations!) using the LinearSolve command. If we write the input as fractions, Mathematica will do the computation exactly and give you the answer as a rational number.

ExactA = {{1036/10000,2122/10000},
               {2081/10000,4247/10000}};
Exactb = {7381/10000,9327/10000};

ExactSolution = LinearSolve[ExactA,Exactb]

Writing an exact expression in approximate form

Mathematica stores everything internally in exact form. You can display these expressions as decimals with the N[ --stuff--, --digits--] command:

ExactSolution

N[ExactSolution,6]

N[ExactSolution,10]

N[ExactSolution,50]

Setting up a solver with limited precision

But let’s be careful here. First, we can write down the general solution to a 2x2 linear equation easily enough:

SymbolicSolution = LinearSolve[{{a,b},{c,d}},{e,f}]

Now if we actually evaluate these fractions, it would look like:

MatrixSolution[{{a_,b_},{c_,d_}},{e_,f_}] := {((d e-b f)/(-b c+a d)),((c e-a f)/(b c-a d))};

MatrixSolution[ExactA,Exactb]

But now we want to be a little bit more careful about precision. Let’s define a function which rounds things to n digits; 
this will be a fixed point representation of the number (because we’re not worrying about making sure the leading 
digit is nonzero).

fx[x_,n_] := N[Round[x*(10^n)]/(10^n),n];

fx[10π,22]

Writing out a computation with roundoff

We now want to evaluate our general solution to the 2x2 matrix problem using fx. This looks awful because we have to 
round after every arithmetic operation:

{((d e-b f)/(-b c+a d)),((c e-a f)/(b c-a d))}

becomes:

MatrixSolution[{{a_,b_},{c_,d_}},{e_,f_},n_] := {fx[((fx[fx[fx[d,n] fx[e,n],n]-fx[fx[b,n] fx[f,n],n],n])/(fx[-fx[fx[b,n] fx[c,n],n]+fx[fx[a,n]fx[ d,n],n],n])),n],fx[((fx[fx[fx[c,n] fx[e,n],n]-fx[fx[a,n]fx[ f,n],n],n])/(fx[fx[fx[b,n] fx[c,n],n]-fx[fx[a,n]fx[ d,n],n],n])),n]};

Solutions with limited precision

We can now check our solutions with limited precision

N[ExactSolution,10]

Now our inputs were specified to four digits of precision, so we expect that we can solve in 4 digits of precision and get a solution that’s accurate to 4 digits

MatrixSolution[ExactA,Exactb,4]

What happened?

One way to understand the problem is to graph approximate solutions to each equation:

f[x_,y_] := 0.1036 x + 0.2122 y - 0.7381
g[x_,y_] := 0.2081 x + 0.4247 y - 0.9327

Here are regions where -(10^(-5))≤ f(x,y) ≤ (10^(-5)) (in yellow) 
and -(10^(-5))≤ g(x,y) ≤ (10^(-5)) (in blue).

Block[{d=4},Show[{RegionPlot[{-(10^(-d)) < f[x,y]< (10^(-d)) ,-(10^(-d)) ≤ g[x,y] < (10^(-d))},
{x,ExactSolution[[1]]-20*(10^(-d)),ExactSolution[[1]]+20*(10^(-d))}, 
{y,ExactSolution[[2]]-20*(10^(-d)),ExactSolution[[2]]+20*(10^(-d))},
PlotStyle->{Yellow,LightBlue}, PerformanceGoal->"Quality",MaxRecursion-> 5,
FrameTicksStyle->Directive[Black,FontSize->16]],
ListPlot[{Callout[ExactSolution,"Exact Solution"]},
LabelStyle->Large]},ImageSize->Large]
]

The lines are almost parallel!

Better demonstration

Here’s a more general demonstration of the solution region for two lines. Consider

x - y  = 0
x - a y = 0

The solution is clearly (0,0) regardless of the value of a. 
But the difficulty of computing the solution numerically
depends on a!

ϵ = 0.001;
Manipulate[RegionPlot[{
-ϵ < x - y < ϵ,
-ϵ < x - a y < ϵ
},{x,-10 ϵ,10ϵ}, {y,-10ϵ,10ϵ},
PerformanceGoal->"Quality",MaxRecursion->5],{a,-1,1}]

Moral of the story

Systems of linear equations such as 

0.1036 x + 0.2122 y = 0.7381
0.2081 x + 0.4247 y = 0.9327

are known as ill-conditioned systems. They are very hard to solve numerically! 
To understand where the algebra goes wrong, look at the symbolic solution:

{((d e-b f)/(-b c+a d)),((c e-a f)/(b c-a d))}

and consider the denominator:

denom =(-b c + a d) /. {a->0.1036,b->0.2122,c->0.2081,d->0.4247}

If we represent it with only 4 digits, we get

rounded = fx[denom,4]

which differs from the right answer by

Abs[(denom-rounded)]/Abs[denom ]

More than 25%!