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

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

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

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

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

The determinant computation is a matter of memory use in terms of how much we want to store for subdeterminants of a Laplace expansion. Mathematica simply refuses to go that route after 11x11. YOu can do your own as below. myDet[mat_] /; Length[mat] <= 4 := Det[mat] myDet[mat_] := myDet[mat] = Sum[mat[[1, j]]*myDet[Drop[mat, {1}, {j}]], {j, ...

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

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 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]]] &, ... 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 This will not answer your question directly, but seems generally on-topic. I have replied to a similar question before on StackOverflow, producing a complete and general Gram-Schmidt process implelentation which I believe to be rather idiomatic for Mathematica. I then specialized it to polynomials. Since it is not going to be migrated, I will simply ... 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 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 First, you need to keep your \thetas and \betas straight. Let's define the matrix in terms of \theta and worry about the relationship with \beta in a bit. epss = -13.6; epsso = -29.1; epsp = -14.1; sssigma = -7.20; spsigma = 9.46; Hmatrix0[theta_] = { {epss, 0, sssigma, Cos[theta]*spsigma, -Sin[theta]*spsigma, 0}, {0, epss, sssigma, ... 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 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 ...

5

If the purpose is as in the following, The reason I want to do that is to filter the list to find the matrix closest to the target one, then an alternative would be to use Nearest: SeedRandom[2013]; matlist = RandomInteger[{0, 4}, {100, 3, 2}]; dist[m1_, m2_] := Norm[m1 - m2, "Frobenius"]; nf = Nearest[matlist, DistanceFunction -> dist]; nf[{{3, ...

5

You can do this using built-in functions using Orthogonalize[x^Range[5], Integrate[#1*#2, {x, -\[Pi], \[Pi]}] &] By hand, a literal implementation of your equations is inn = Integrate[#1*#2, {x, -\[Pi], \[Pi]}] &; v = x^Range[5]; u = ConstantArray[0, Length@v]; u[[1]] = v[[1]]; Do[ u[[i]] = v[[i]] - Total@Table[ u[[j]]*inn[v[[i]], ...

5

Prior to version 5 users had to load an add on package that contained a tridiagonal solver based on the Thomas algorithm I think the code for those old packages is probably accessible somewhere. I found this in the library archive: http://library.wolfram.com/infocenter/MathSource/4827/ This package uses Compile. The earlier built in add-on package was ...

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

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

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] + ...

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