# Tag Info

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 = ...

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

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

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]:= ...

11

Looking at CompilePrint[compiledGlynnAlgorithm] there are some CopyTensor in it which aren't really needed. There's also a few CoerceTensor in there when it might be faster to just coerce the integer matrix once at the beginning. By slightly adjusting the function all CopyTensor and CoerceTensor go away giving a small increase in speed: ...

11

I had this laying around from a course in numerical linear algebra I taught a few years ago. Here's a matrix whose nonzero elements describe the basic shape. size = 50; nw = Partition[Table[i, {i, 1, size^2}], size]; sw = Partition[Table[i, {i, size^2 + 1, 2*size^2}], size]; se = Partition[Table[i, {i, 2*size^2 + 1, 3*size^2}], size]; L = ...

11

You can always modify the matrix so that the most negative eigenvalue is also the one with the largest absolute value, and hence corresponds to the first in the list returned by EigenVectors. An upper bound for the largest absolute value of any eigenvalue is the Hilbert-Schmidt norm. So you can rescale your matrix by subtracting this norm times the unit ...

10

Working with LinearSolve we encounter some inconsistency of the related option Modulus -> z if z is not prime. Nonetheless we could do this Mod[ LinearSolve[ {{1, 1, 1}, {4, 2, 1}, {9, 3, 1}}, {31, 3, 11}], 54] {18, 26, 41} Unfortunately we can get only one solution unlike when working with Solve. These posts describe another problems or bugs ...

9

As $P$ is explicitly constructed from eigenvectors of a self-adjoint matrix, it is unitary, i.e $P P^\dagger = I\qquad$ where the $\dagger$ is the conjugate transpose (or Hermitian conjugate, if you prefer). So, calculating the inverse is simply ConjugateTranspose[P] which is much faster than calculating it using Inverse. That said, you have to ensure that ...

9

I'm not sure one can get the larger cases due to memory needs. The code below, which uses the array generation of @rcollyer, seems reasonably effective at least for the smaller end. n = 2500; Timing[sum = 0; diags = Sqrt[ 1./HarmonicNumber[n, 2]]/(Sqrt[n - #] (# + 1) &@Range[0., n - 1]); hh = Array[If[#2 < #1, 0., diags[[#2 - #1 + 1]]] &, ...

9

Here is a way to get the eigenvectors, using NullSpace. I'm using a smaller matrix as an example here: MatrixForm[matrix = N[HilbertMatrix[4]]] $$\left( \begin{array}{cccc} 1. & 0.5 & 0.333333 & 0.25 \\ 0.5 & 0.333333 & 0.25 & 0.2 \\ 0.333333 & 0.25 & 0.2 & 0.166667 \\ 0.25 & 0.2 & 0.166667 & 0.142857 ... 9 From the Documentation on Numerical Implementation: Machine-precision matrices are typically converted to a special internal representation for processing. SparseArray with rules involving patterns uses cylindrical algebraic decomposition to find connected array components. Sparse arrays are stored internally using compressed sparse row formats, ... 7 I verified that Mathematica returns the correct set of eigenvalues and eigenvectors for this matrix by comparing them to MATLAB's output. I'll show how to do the comparison, as this might help reveal mistakes in your own comparison, if there were any. (Note that if there are degenerate eigenvalues, then the eigenvectors are not unique, so there may be ... 7 64 bit Mathematica does not have any practical limits on this. What limits you is the speed of your computer and the available memory. A k\times k matrix will take a bit more than 8\times k^2 / 1024^3 gigabytes of memory, so you see that a 10^6 \times 10^6 matrix needs ~7500 GB of memory to store. You probably don't have that much in your computer. ... 7 The determinant of the matrix A in this case is about 10^282 The determinant isn't very useful, but the condition number is: You can use SingularValuesList to get the largest and the smallest singular value. If the ratio between the two is too large, the matrix is ill-conditioned. Solving an ill-conditioned linear system will still give "exact" results ... 7 You might get a speed up by restricting compiledGlynnAlgorithm to work on just one row of the Gray Code list, allowing the Listable and Parallelization to come into play. I say "might" because the speed up will depend on the details of your hardware. Redefine compiledGlynnAlgorithm like so (note that it now takes a one dimensional list for d): ... 6 The idea behind the Jordan normal form does the trick, even though JordanDecomposition does not. (Incidentally, this suggests there may be a more reliable, stable algorithm to obtain Jordan decompositions than is implemented in Mathematica...) The resulting solution is very short, efficient, and numerically stable when applied to floating-point matrices. ... 6 StateSpaceModel will linearize the equations. The system is nonlinear. So let's look at the nonlinear solution first. The parameters: a1 = {{-3, 2}, {-0.25, 1}}; a2 = {{-1.9, -0.4}, {-2.24, -4.7}}; b1 = {{0.25}, {1}}; b2 = {{-2.5}, {1}}; c = {{1, 0.5}, {0, 1}}; Subscript[ρ, 1] = (1 - Tanh[Subscript[x, 1][t]])/2; Subscript[ρ, 2] = 1 - Subscript[ρ, 1]; Set ... 6 Higher-order SVD (in sense of Tucker decomposition) of the matrix M with dimensions d_1\times d_2\times\cdots\times d_n is$$ M_{i_1,i_2,\dots,i_N} = \sum_{j_1} \sum_{j_2}\cdots \sum_{j_N} s_{j_1,j_2,\dots,j_N} u^{(1)}_{i_1,j_1} u^{(2)}_{i_2,j_2} \dots u^{(N)}_{i_N,j_N},  where $s$ is the core tensor and $u^{(i)}$ is the orthogonal matrix. The matrix ...

6

The symbol $\epsilon$ is called the Levi-Civita tensor, it is usually written in terms of its components $\epsilon_{i j k}$. Having defined two vectors (i.e. first order contravariant tensors): a = { a1, a2, a3}; b = { b1, b2, b3}; we have to build up a tensor product of two vectors and Levi-Civita tensor, i.e. we get a new tensor of valence $(3,2)$ ...

5

Update I got a MatrixRank of 4 with the original approximate data, but with the updated exact data, the rank is 3. The basic idea is that Orthogonalize will return an orthonormal basis for the subspace spanned by the vectors, along with some zero vectors interspersed. (Orthonormal means unit length vectors that are pairwise perpendicular.) Deleting the ...

5

Well, you certainly can if you have version 9 and don't mind using the (quite verbose) tensor notation: Assuming[ {b ∈ Complexes, (x | y) ∈ Matrices[{n, n}]}, TensorReduce@ TensorContract[ b TensorContract[TensorProduct[x, y, y], {{2, 3}, {4, 5}}] + b TensorContract[TensorProduct[y, x, y], {{2, 3}, {4, 5}}], {{1, 2}} ] ] (* -> 2 b ...

5

Yes, you can simplify the trace this way. First, Tr[x.y.z] is invariant under cyclic permutations, so Tr[x.y.z] = Tr[y.z.x] = Tr[z.x.y], as described in Wikipedia's entry on trace of a product. For the case of interest, this means Tr[y.x.y] is equal to Tr[x.y.y]. Hence Tr[b x.y.y + b y.x.y] = Tr[b x.y.y] + Tr[b y.x.y] = b Tr[x.y.y] + ...

5

Not too hard; all that's needed is a simple application of matrix identities: ColumnHermiteDecomposition[mat_ /; MatrixQ[mat, IntegerQ]] := Transpose /@ HermiteDecomposition[Transpose[mat]] Test: mat = {{1, 2, 3, 2, 2}, {1, 2, 3, 4, 0}, {0, 5, 4, 2, 1}, {3, 2, 4, 0, 2}}; {u, t} = ColumnHermiteDecomposition[mat]; u {{8, 24, 22, ...

5

We can easily verify the assumption that your matrices are equal is wrong, e.g. let's take two 2 x 2 matrices: A = {{10, 5}, {2, 2}}; B = {{1, 0}, {1, 6}}; now we have A.B.Transpose[A] == A.Transpose[B].Transpose[A] False However a simple fact in linear algebra says that the former matrix is equal to the latter one transposed. TraditionalForm[ ...

5

Update: re-run the tests for 10 times each, to get better average readings. (and also corrected a coding error) Conclusion: Default and \$MachinePrecision are fastest. (hard to see any difference, did 10 times per each). So if you want fast, just use Default. This page talks more about eigenvalue computation in Mathematica and the use of Lapack and under ...

5

By default the computation will be done in machine precision, without precision tracking. I believe this is the fastest method you can get without some form of packing that places multiple values in a single machine float, which I know little about. Once you use SetPrecision you are engaging the arbitrary precision engine with precision tracking, which is ...

5

If the vectors v come one at a time and must be handled that way, as opposed to being collected for batch analysis, then I don't see how you can beat simply v.B.v. However, if you can collect them as the rows of a many x n matrix V then Total[V.B*V,{2}] is very fast. (If the vectors come "naturally" as the columns of an n x many matrix U then use ...

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