# Tag Info

32

Mathematica is every bit as fast as Matlab for these types of computations. The source of the discrepancy arises from the fact that Timing keeps track of total time used by all processors when Mathematica distributes the computation across them. We can examine a fair comparison using AbsoluteTiming, which is more comparable to Matlab's tic and toc. ...

23

You're looking for ArrayFlatten. For your example matrices, R = ArrayFlatten[ {{A, {t}\[Transpose]},{0, 1}} ] (* => {{1, 0, 0, 1}, {0, 0, 1, 1}, {0, -1, 0, 1}, {0, 0, 0, 1}} *) The construct {t}\[Transpose] is necessary for ArrayFlatten to treat t as a column matrix. Then to find $\boldsymbol{R}^{-1}$, you run Inverse[R] (* => {{1, 0, 0, ...

20

(I've been waiting for somebody to ask this question for months... :D ) Here's the Mathematica implementation of the Frobenius companion matrix approach discussed by Jim Wilkinson in his venerable book (for completeness and complete analogy with built-in functions, I provide these three): PolynomialEigenvalues[matCof : {__?MatrixQ}] := Module[{p = ...

17

Mathematica does not support this directly. You can do things of this sort using an external package called NCAlgebra. http://math.ucsd.edu/~ncalg/ The relevant documentation may be found at http://math.ucsd.edu/~ncalg/DOWNLOAD2010/DOCUMENTATION/html/NCBIGDOCch4.html#x8-510004.4 In particular have a look at "4.4.8 NCLDUDecomposition[aMatrix, Options]" ...

17

This answer is almost entirely about mathematics and algorithms, not Mathematica implementation. I'm not sure whether such answers are welcome on this site; I hope I haven't offended. This is a very important problem in computational algebra, but is usually stated in a more sophisticated way. The usual way that one thinks about it is to consider the ...

16

It's just a matter of the difficulty inherent in the numerical computation of determinants. Here's what Cleve Moler has to say about determinants and characteristic polynomials in chapter 10 of his book on numerical computing: Like the determinant itself, the characteristic polynomial is useful in theoretical considerations and hand calculations, but ...

15

acl already posted the crucial information needed to solve this conundrum (i.e., the definition of InternalCompileValues[LinearSolve]), but wishes to delete his post since he had not interpreted it to give the complete answer. Therefore I re-post the following observation along with a summary of what it means. The input, InternalCompileValues[]; ...

15

Here is the method I outlined. I'll illustrate on a small example where we split matrix into top and bottom halves. In[794]:= SeedRandom[1111]; halfsize = 3; mat = RandomInteger[{-4, 4}, {2*halfsize, 10}] Out[796]= {{-3, -1, 3, -3, 3, 3, 3, 3, 4, 2}, {3, 3, -3, 0, 0, 1, -2, -4, 0, -1}, {-3, 4, 3, 0, -2, 4, 3, -2, -2, -2}, {2, 2, 4, 0, -4, 4, -1, -4, ...

15

Mathematica offers a pretty complete set of functionality for linear algebra, and it has improved in recent versions. For example, since version 5, Mathematica has offered the generalised Schur decomposition (also known as the QZ decomposition). This certainly wasn't available in earlier versions. It handles sparse matrices and many other wrinkles. And if ...

15

Here's my attempt. To get the matrix representing the Laplacian I use LaplacianFilter on an array of symbols and CoefficientArrays to extract the coefficients. n = 200; shape = ArrayPad[ConstantArray[0, {n/2, n/2}], {{0, n/2}, {0, n/2}}, 1]; shapeVector = Flatten @ Position[Flatten @ shape, 1]; symbolArray = Array[x, {n, n}]; symbolLaplacian = ...

14

Based upon your update, you are trying to solve the system $$\mathbf{A}\vec{x} = \vec{b}$$ for $\vec{x}$, so LinearSolve is exactly what you want. Also, it has the exact form LinearSolve[A, b] that you're asking for. Internally it uses a form of Gaussian elimination to solve such systems; this is most likely a variant of LU decomposition, but other ...

14

Actually what you want is HermiteDecomposition. It is the integer ring form of RowReduce; that latter, while working largely over the integers for the forward elimination, actually is an echelon form over the rational field. As for efficiency, you'll have to experiment to see if it meets your needs. If not, feel free to post or send examples so I'll have ...

14

The rank of a matrix is typically determined by performing a Gaussian elimination and is given by the number of non-zero rows. In your second case, the large number $5.4\times 10^{12}$, when eventually used as a pivot, gives a badly conditioned matrix (myM2 is the second matrix in your question): RowReduce[myM2] RowReduce::luc: Result for RowReduce of ...

13

We can construct this matrix directly as a SparseArray. This allows some classes of numerical matrices to be stored as packed arrays while being combined with symbolic or exact vectors (or vice versa), so there can be storage and run-time efficiency reasons for using a SparseArray, in addition to the obvious benefit of direct construction. On the other hand, ...

13

Here is a very simple way to do it: Table[1/i! D[M, {a, i}] /. a -> 0, {i, 0, 3}] (* ==> {{{15, 0}, {0, 2}}, {{0, 1}, {1, 0}}, {{1, 5}, {-5, 0}}, {{0, 0}, {0, 0}}} *) This works even if the entries are not polynomials. If they are, you can replace the arbitrary maximum 3 in the Table index by the degree of the polynomial: Max[Exponent[M, a]] ...

13

Use Eigenvalues[mat, Cubics -> True] Eigenvectors[mat, Cubics -> True] sometimes Quartics -> Truecan be needed. or ToRadicals @ Eigenvalues[ mat] ToRadicals @ Eigenvectors[ mat] In general one cannot find roots (of higher order) polynomials in terms of radicals. The reason that Mathematica allows this option is that in general it is ...

13

This might be as good a time as any to distill the collective wisdom of Messrs. Huber, McClure, and Toad R. M. As already mentioned, there is this quantity of great interest to people in the business of solving simultaneous linear equations, called the condition number, and conventionally denoted by the symbol $\kappa$. This is usually associated with a ...

13

There is no need to use Eigensystem or Eigenvectors to find the axis of a rotation matrix. Instead, you can read the axis vector components off directly from the skew-symmetric matrix $$a \equiv R^T-R$$ In three dimensions (which is assumed in the question), applying this matrix to a vector is equivalent to applying a cross product with a vector made up ...

13

There is a lot going on with your code, as already pointed out by others. But, I would like to point out a few more things, and this was much, much to long for a comment. First, I would use SparseArray to generate H, instead of the pre-allocate, fill-in method you are using, as follows: H = SparseArray[ {l_, m_} /; m >= l :> Sqrt[pp]/(Sqrt[n + l - ...

12

The keyboard commands Ctrl+Enter, Ctrl+, and Tab can be used to enter this format. You can also use the menu Insert > Table/Matrix to create a table of specified size with placeholders. See Entering Tables and Matrices. Depending on the meaning of the question, this may have some bearing:

12

Use the following: ListConvolve[a, b, {1, -1}, 0] (* ==> {1, 3, 6, 10, 10, 10, 9, 7, 4} *) This says: align the first element of b with the 1st element of a, align the last element of b with the -1st (i.e. last) element of a, pad with 0s if necessary. Have you tried searching the docs for "convolution"?

12

By default, Mathematica assumes symbols to be complex. However the elements on the main diagonal of a Hermitian matrix are necessarily real. To force Mathematica to interpret the elements on diagonal of m to be real you could replace them by their real part, i.e. m = {{Re[n], a, b, b}, {Conjugate[a], Re[n], b, b}, {Conjugate[b], Conjugate[b], Re[c], ...

12

Here is simple (unweighted) Mma version of the Matlab implementation of Covariance Bending. ClearAll[covBending]; covBending[mat_, tol_: 1/10000]:=If[PositiveDefiniteMatrixQ[mat], mat, NestWhile[(Eigensystem[#][[2]].DiagonalMatrix[ Max[#, tol] & /@(Eigensystem[#][[1]])].Transpose[ Eigensystem[#][[2]]]) &, N@mat, Min[Eigensystem[#][[1]]] < ...

12

Using Replace (assuming you only want to replace on level 2, as you mention "matrix"): Using Except My first version (I kept this version to point out the usage/impact of Orderless) Replace[SIGMA, Except[HoldPattern[___ x^2 ___]] -> 0, {2}] {{0, b x^2}, {0, 0}} Improved version, thanks to Leonid Replace[SIGMA, Except[___ x^2] -> 0, {2}] as ...

12

While the other answers are nice, the icon deserves a closer look: Note, in particular, that four of the six edges are not constrained by the ostensible Dirichlet boundary conditions, nor is it clear that they solve a Neumann problem. And indeed, as I noted in the comments this is supported by the OP's first link. In short, to produce the logo, they took ...

11

Initially, Mathematica is not designed for such abstract calculations. But, Mathematica is a powerful programming language, so that one can add such functionality easily. See the following examples in related area of differential geometry: calculations in symbolic dimensions Abstract calculations

11

If you need to be sure that the order is correct, there is a function Eigensystem that returns a list of both the eigenvalues and -vectors, which is in the right order. {eValues, eVectors} = Eigensystem[{{2, 0}, {0, 1}}]; eValues eVectors {2, 1} {{1, 0}, {0, 1}} It's probably worth using just for the slight off-chance that Eigenvalues and Eigenvectors ...

11

A reasonable alternative, is to construct an explicitly Hermitian matrix by exploiting the fact that any matrix, $M$, can be written as the sum of a Hermitian matrix, $H$, and a skew-Hermitian (or anti-Hermitian, if your in physics) matrix, $S$. This implies that a matrix can be made Hermitian simply by $$H = \frac{1}{2}(M + M^\dagger)$$ or skew-Hermitian ...

11

Since you're working with vectors, just let Mathematica know that these are vectors. Some other systems (MATLAB and its relatives in particular) have the limitation that they can only work with matrices, forcing you to distinguish between row vector and column vectors and keep transposing. This is not necessary nor convenient in Mathematica. In[1]:= ...

Only top voted, non community-wiki answers of a minimum length are eligible