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


15

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


13

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


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


12

With the understanding that the criterion for including a point in the 3D matrix is that it has 3 unique domains in any of the 26 points surrounding a zero value in the watershed here is a simple way to extract the data. In the case where zero values in the watershed are in minority it may be the fastest approach to getting a list of points that fulfill the ...


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

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


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


8

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


7

You have to be a bit careful here; your last approach does not give the desired matrix F: (LinearSolve[dV, dv] // Transpose).dV === dv False Then what does give the correct output? We can use the fact that for generic matrices $A$ and $B$ we have $A.B = (B^T.A^T)^T$ and write F = Transpose @ LinearSolve[Transpose @ dV, Transpose @ dv] And indeed, ...


6

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


6

I'll make this a response since it seems to be a popular topic today (it showed up independently on the Wolfram Community forum). One can usually do this by finding the eigenvectors for a random linear combination of the matrices. This works for nonderogatory matrices, that is, ones that do not have nontrivial blocks in their Jordan decomposition. In ...


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 compatibility information at Compatibility/tutorial/LinearAlgebra/MatrixManipulation says These functions were available in previous versions of Mathematica and are now available on the web at library.wolfram.com/infocenter/MathSource/6770: LinearEquationsToMatrices InverseMatrixNorm ConditionNumber You can download the original package ...


6

The order of eigenvalues is the most convenient order for the algorithm, which find these eigenvalues. You can always order them as you want very simply a = # + #\[Transpose] &@RandomReal[1, {10, 10}]; {ε, ψ} = Eigensystem[a]; {ε, ψ} = {ε[[#]], ψ[[#]]} &@ Ordering[ε]; Furthermore, the eigenvalues can be complex for non-Hermitian matrices. There is ...


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


5

The set of equations given by $$e_p = \sum_{i=1}^D c_i\ p^i,\quad p\in\{1,\ldots,D\}$$ where the $e_p$ are the diagonal elements (or eigenvalues) of $\mathbf{E}$ can be written as $$\mathbf{e} = \mathbf{c}\ \left( \begin{array}{ccc} 1 & 1 & \ldots & 1^{D-1} \\ 1 & 2 & \ldots & 2^{D-1} \\ \vdots & & \ddots & \vdots ...


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

You could define a function that constructs the product of Pauli matrices as follows. I use KroneckerProduct here because you are planning to form the matrix product with a $2L\times2L$ matrix, so we have to have the Pauli matrices arranged in a corresponding block matrix: pauliProduct[n_] := Module[{l = Length[n]}, Total@MapIndexed[ KroneckerProduct[ ...


5

The two eigenvalues are degenerate. In general whenever you have degenerate eigenvalues there is arbitrariness in how to select the eigenvectors. You can construct a constant vector from the two degenerate eigenvectors like so: m = Import["https://dl.dropboxusercontent.com/u/63413473/ExampleMatrix.mx"]; {evals, evecs} = Eigensystem[m, -2]; You can check ...


5

The problem here can arise because of numerical underflow which appears for sufficiently large dimension of the problem. Some numerically very small number multiplies the parameter "a" and therefore "a" does not appear in the "solution". Consider a simple example Define the matrix m (fill it with random numbers, here exponentially distributed) In[263]:= ...


5

Dense blocks One can split matrices by blocks and use these blocks as a dense matrices blockSize = 100; m1[a_] m2[b_] ^:= a.b; part = Developer`PartitionMap[If[Length@#@"NonzeroValues" > 0, m1@Normal@#, 0] &, a, {blockSize, blockSize}]; a2.u_ ^:= Flatten[Developer`ToPackedArray[ part.Developer`PartitionMap[m2, u, blockSize]], 1] I have ...


5

Ok, I'll go ahead and answer my own question. After some thinking I came up with this: getIntersection[l1_, l2_] := Module[{n, ker, coeffs}, ker = NullSpace[Transpose[Join[l1,l2]]]; n = Length[l1]; coeffs = Map[Function[v, v[[1 ;; n]]], ker]; Return [Map[Function[v, v.l1], coeffs]]; ]; Can this be made any faster?


5

Maybe MapThread gives the levelspec control you need: lst = {{{a, b}, {c, d}}, {{w, x}, {y, z}}}; g @@ MapThread[f, lst, 1] (* or just g @@ MapThread[f, lst] -- thanks: Mr.W *) (* g[f[{a,b},{w,x}],f[{c,d},{y,z}]] *) g @@ MapThread[f, lst, 2] (* g[{f[a,w],f[b,x]},{f[c,y],f[d,z]}] *) or, g @@ Thread[f @@ lst] (* g[f[{a,b},{w,x}],f[{c,d},{y,z}]] *) or ...


4

Suppose we work with objects in some symbolic dimension dim: In[1]:= $Assumptions = (e[1] | e[2] | e[3] | e[4]) \[Element] Vectors[dim] && (a | b) \[Element] Matrices[{dim, dim}, Antisymmetric[{1, 2}]]; Now we can do something like In[2]:= TensorRank[e[1]\[TensorWedge]e[2]] Out[2]= 2 and basic manipulations like In[3]:= ...


4

note see bottom of answer for a required fix for version 10+ A 2D example of what i suggested in my comment: Based on Belisarius example..we get three lines obviously not connected properly: m = SparseArray[{i_, j_} -> Sin[i j 9 /10 y], {3, 3}]; alle = Table[Eigenvalues[m], {y, 0, 1, .1}]; original = Show[ MapIndexed[ListPlot[Flatten[Take[alle, ...


4

Basically it amounts to using CoefficientArrays along with a bit of massaging for sense of inequality. getLinearForm[ineqs : {(_LessEqual | _GreaterEqual) ..}, vars_] := Module[ {mults, lpolys, rhs, lhs}, mults = Map[Head, ineqs] /. {LessEqual -> 1, GreaterEqual -> -1}; lpolys = mults*Apply[Subtract, ineqs, {1}]; {rhs, lhs} = ...



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