Group Abstract Group Abstract

Message Boards Message Boards

0
|
8.3K Views
|
8 Replies
|
4 Total Likes
View groups...
Share
Share this post:

Numerical Error with Matrix operations

Gordon
Posted 11 years ago
POSTED BY: Gordon
8 Replies
Sort By:

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 BY: Guillermo Hernández
Gordon
Gordon
Posted 11 years ago

It's so good to learn ArrayFlatten. Helps a lot.

Thank you
POSTED BY: Gordon
POSTED BY: Daniel Lichtblau
Gordon
Gordon
Posted 11 years ago

Thanks for pointing out my mistake on [A|b] construction. That's the best answer.
POSTED BY: Gordon

My apologies, I forgot about the last column as the title mentioned 'numerical errors'. Thanks for completing the answer.
POSTED BY: Guillermo Hernández

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 BY: Ivan Morozov
Gordon
Gordon
Posted 11 years ago

Wow, it works. Why? What's the magic behind that?
POSTED BY: Gordon

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.
POSTED BY: Guillermo Hernández
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard