0
|
6762 Views
|
8 Replies
|
4 Total Likes
View groups...
Share
GROUPS:

# Numerical Error with Matrix operations

Posted 10 years ago
 A is a 3x3 matrix, b is a 3x1 vector. I try to convert [A|b], a 3x4 matrix, to [I|0]. So the formula is right multiple [A|b] by a matrix: $$\begin{bmatrix} A^{-1} & -A^{-1}b\ 0 & 1\ \end{bmatrix}$$ Here's the code in Mathematica 10 A=Table[RandomReal[],{i,3},{j,3}] b=Table[RandomReal[],{i,3}] M=Transpose[Join[A, {b}]] ret=M.Join[Transpose[Join[Inverse[A], {-Inverse[A].b}]], {{0, 0, 0, 1}}] The last column of ret is always non-zero. The ret will often look like for example $$\begin{bmatrix} 1. & 0. & 0. & -0.979237\ 0. & 1. & 4.44089*10^{-16} & 2.17005\ -8.8817810^{-16} & -5.5511210^{-17} & 1. & -1.79119 \end{bmatrix}$$ How can I make Mathematica do it right?
8 Replies
Sort By:
Posted 10 years ago
 A straightfoward method to check the matrix product you asked for: define your matrix by blocks using ArrayFlatten: A = Table[RandomReal[], {i, 3}, {j, 3}]; b = List /@ Table[RandomReal[], {i, 3}]; (*So it is a matrix*) M = ArrayFlatten[{{A, b}}]; ret = ArrayFlatten[{{Inverse[A], -Inverse[A].b}, {0, 1}}]; (*Notice numbers represents matrices for ArrayFlatten*) M.ret // Chop // MatrixForm Which returns the 3x4 identity matrix we were expecting.
Posted 10 years ago
 It's so good to learn ArrayFlatten. Helps a lot.Thank you
Posted 10 years ago
 I guess I should state that Mathematica is doing it right. Truncation error is not avoidable with finite precision arithmetic. As for the problem at hand, there are a couple of things that seem to be amiss.(1) Create the augmented matrix {A|b] correctly. This can be done with M = Join[A, List /@ b, 2]. Your version transposes A in the process.(2) Now just use RowReduce on it. No inverses, no multiplications.Here is an example. SeedRandom[1111]; aA = Table[RandomReal[], {i, 3}, {j, 3}] b = Table[RandomReal[], {i, 3}] mat = Join[aA, List /@ b, 2] rred = RowReduce[mat] Out[434]= {{0.0777960697813, 0.556265936296, 0.807681686834}, {0.527982841158, 0.555876762387, 0.981090055861}, {0.907596678267, 0.364390115148, 0.722450567448}} Out[435]= {0.420186553434, 0.1560084863, 0.13002033209} Out[436]= {{0.0777960697813, 0.556265936296, 0.807681686834, 0.420186553434}, {0.527982841158, 0.555876762387, 0.981090055861, 0.1560084863}, {0.907596678267, 0.364390115148, 0.722450567448, 0.13002033209}} Out[437]= {{1, 0., 0., 0.415692234671}, {0, 1, 0., 4.46163163168}, {0, 0, 1, -2.59261340179}}Check that this agrees with LinearSolve (I assume that was the purpose). rred[[All, -1]] == LinearSolve[aA, b] Out[438]= True 
Posted 10 years ago
 Thanks for pointing out my mistake on [A|b] construction. That's the best answer.
Posted 10 years ago
 My apologies, I forgot about the last column as the title mentioned 'numerical errors'. Thanks for completing the answer.
Posted 10 years ago
 Gordon, Replace -Inverse[A].b with -b.Inverse[A] and then use Chop as suggested by Guillermo. A = Table[RandomReal[], {i, 3}, {j, 3}] ; b = Table[RandomReal[], {i, 3}] ; M = Transpose[Join[A, {b}]] ; ret = M.Join[ Transpose[Join[Inverse[A], {-b.Inverse[A]}]], {{0, 0, 0, 1}}] ; ret // MatrixForm ret // Chop // MatrixForm I.M.
Posted 10 years ago
 Wow, it works. Why? What's the magic behind that?
Posted 10 years ago
 Remember that when you perform numerical computations with floating point numbers they have a limited precision, which might lead to numerical errors. It is usual to assume in such cases that results that are very close to zero are zero indeed. Mathematica's Chop function can be used for that. Just try: Chop[ret] You can also check it is a precision error by using arbitray precision for your numers: A=Table[RandomReal[WorkingPrecission->30],{i,3},{j,3}] b=Table[RandomReal[WorkingPrecission->30],{i,3}] M=Transpose[Join[A, {b}]] ret=M.Join[Transpose[Join[Inverse[A], {-Inverse[A].b}]], {{0, 0, 0, 1}}] Which will make those almost-zeros even smaller.You can also check this tutorial for more information.