Tag Info

Hot answers tagged

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


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


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

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


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

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

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

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

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): ...


7

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)$ ...


7

s = {x, y} /. Solve[a x + y == 7 && b x - y == 1, {x, y}][[1]] {8/(a + b), -((a - 7 b)/(a + b))} lsa = LinearSolve[{{a, 1}, {b, -1}}, {7, 1}] {8/(a + b), (-a + 7 b)/(a + b)} f = LinearSolve[{{a, 1}, {b, -1}}]; lsb = f[{7, 1}] // Simplify {8/(a + b), -((a - 7 b)/(a + b))} s == lsa == lsb // Simplify True Solve can handle a ...


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

Unfortunately the various Solve like functions don't show the same behavior concerning string-type variable names. Some do accept strings (e.g. NDSolve in version 9): NDSolve[{"x"'["t"] == 0.1*"x"["t"], "x"[0] == 1}, "x", {"t", 0, 1}] but others don't. For your case, Eliminate obviously doesn't, but you only need to convert those variables which you want ...


6

x = {3, 1} y = {2, 5} a Defer[x] + b Defer[y] /. First@Solve[a x + b y == {7, 11}, {a, b}] (* x + 2 y *) Note that the output is usable as such: evaluate it and you'll get the combination result. If only integer coefficients are desired, changing the Solve to something like: Solve[a x + b y == {7, 11} && {a, b} \[Element] Integers, {a, b}] ...


6

The phase (and length) of the eigenvectors is completely undetermined unless you specify extra conditions in addition to the eigenvalue equation. Given that you don't have any additional conditions, it's not surprising that there is no well-defined way to plot the real and imaginary parts of each eigenvector component. A simple condition that makes the ...


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

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

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

Something along these lines perhaps? coefs[pol_, base_] := Last@CoefficientArrays[pol, base] Star[lhs_, rhs_] := coefs[lhs, starBase].Transpose /@ starTen.coefs[rhs, starBase].starBase So you could define starBase = v~Array~4; starTen = Array[c, {4, 4, 4}]; and now Star[5 v[3] + 7 v[4], 14 v[1]] (* 14 (5 c[3, 1, 1] + 7 c[4, 1, 1]) v[1] + ...


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


5

I don't think you need to use any specific tensor functionality. SolveAlways seems to suffice: SolveAlways[ T1 - T3 == a1 T1 + a2 T2 + a3 T3, {T1, T1, T3} ] (* => {{a1 -> 1, a2 -> 0, a3 -> -1}} *)


5

Given vectors $\;v_1, v_2, Tv_1, Tv_2$: { v1, v2} = {{2, -1}, {1, -1}}; {Tv1, Tv2} = {{4, -1}, {5, -3}}; This is a direct way of solving underlying linear system: With[{ m = Array[a, Dimensions[{v1, Tv1}]]}, m.Subscript @@@ {{x, 1}, {x, 2}} /. Solve[ m.#1 == #2 & @@@ {{v1, Tv1}, {v2, Tv2}}, Join @@ m] // Transpose // ...


5

You can do it with FindGeometricTransform[]: data = {{13.2, 200, 58, 21.2}, {10, 263, 48, 44.5}, {8.1, 294, 80, 31}, {8.8, 190, 50, 19.5}, {9, 276, 91, 40.6}, {7.9, 204, 78, 38.7}, {3.3, 110, 77, 11.1}, {5.9, 238, 72, 15.8}, {15.4, 335, 80, 31.9}, {17.4, 211, 60, 25.8}}; {error, f} = FindGeometricTransform[PrincipalComponents[data], data]; ...


5

Here is a variant adapted from this MathGroup thread permanentC = Compile[{{m, _Real, 2}}, With[{len = Length[m]}, (-1)^len*Module[ {s = {0.}, u = 0.}, Do[ s = N[IntegerDigits[n, 2, len]]; u += (-1)^Round[Total[s]]*(Times @@ (m.s)), {n, 2^len - 1}]; u]], CompilationTarget -> "C"]; I checked it on the test set ...


5

There is an alternative solution to implement AMD into Mathematica, by using the MinCut function from the GraphUtilities Package. This function is not just working on Graphs but is also usable for Matrices. The algorithm is reordering the Matrix into blocks and effectively reducing the off diagonal elements. It can be involved rather easily. If one starts ...



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