Example 2. 

Effect of precision on the solution 
to a system of linear 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

but to make it easier to manipulate in Mathematica, we will put the constants on the left hand side.

To solve the system, we eliminate x by multiplying the first equation by alpha = (0.2081/0.1036) and subtracting it from the second equation.  Mathematica will do this as rational numbers, with no loss of precision. Here is the correct solution:

Eq1F =  (1036/10000)x + (2122/10000) y - (7381/10000)
Eq2F  =  (2081/10000) x + (4247/10000) y - (9327/10000)

alphaF =(Coefficient[Eq2F,x])/(Coefficient[Eq1F,x])

Eq3F= Coefficient[Eq2F,y]*y  - (9327/10000)-(alphaF*Coefficient[Eq1F,y]*y - alphaF*(7391/10000))   

 TrueAnswer = y /. Solve[Eq3F == 0,y]
 TrueXAnswer = N[x /. Solve[ Eq1F ==0 /. Solve[Eq3F == 0,y], x]]

Now we compare the true answer (rounded to MyPrecision digits) to the result obtained by doing all the arithmetic in the solution in MyPrecision digits:

MyPrecision = 3;
p[x_] := (  (Round[(10^MyPrecision)*x]/10^MyPrecision) )

Eq1 =  p[0.1036] x + p[0.2122] y - p[0.7381]
Eq2  = p[0.2081] x + p[0.4247] y - p[0.9327]

alpha =p[(p[Coefficient[Eq2,x]]/p[Coefficient[Eq1,x]])]

Eq3 = p[Coefficient[Eq2,y]]*y  - p[0.9327]-(p[alpha*Coefficient[Eq1,y]]*y - p[alpha*(p[0.7391])])      

 N[p[y /. Solve[Eq3 == 0,y]],MyPrecision]
 N[TrueAnswer,MyPrecision]

What happened? One way to see the problem is to graph the two lines and see what the intersection looks like. Plotting neighborhoods of the lines shows us that the lines are nearly parallel!

Show[{RegionPlot[Eq1F < 0.1 && Eq1F > -0.1,{x,-700,-750}, {y,355,360},PlotStyle->Yellow, PerformanceGoal->"Quality",MaxRecursion-> 5],
RegionPlot[Eq2F < 0.1 && Eq2F > -0.1,{x,-700,-750}, {y,355,360},PerformanceGoal->"Quality", MaxRecursion-> 5]}]

It's actually pretty hard to display the region where both of these equations are approximately satisfied to a smaller tolerance in Mathematica, but here's a demo which shows you how the approximate solutions of two "model" linear equations: x - y = 0 and x - a y = 0 vary as we change the parameter a.

Model1 = x - y;
Model2 = x - a*y;
Epsilon = 0.01;
Manipulate[RegionPlot[{Model1 < Epsilon && Model1 > -Epsilon,x - a y < Epsilon && x - a y > -Epsilon},{x,-10*Epsilon,10*Epsilon}, {y,-10*Epsilon,10*Epsilon}, PerformanceGoal->"Quality",MaxRecursion-> 5],{a,-1,1}]

We see that the problem is the shape of the region of intersection between the approximate solutions of x - y = 0 (blue) and the approximate solutions of x - ay = 0 (pink). 

If a is close to 1, then the system becomes very difficult to solve accurately (even though we know the true solution is (0,0)).