# Tag Info

0

As per comment by @bill, I have not been able to reproduce the problem. I did a trial run using the original matrix using the code below, and the figure values seem fixed. However there are a no. of eigenvectors with low probability of existence (see first figure), and might give rise to problems. ListPlot[Eigenvectors[N[mat]], Frame -> True, ...

3

This is a somewhat high-brow way of showing the Cayley-Hamilton theorem, through the power of holomorphic functional calculus. As I mentioned in this answer, one of the standard ways to define a matrix function is through a Cauchy-like construction: $$f(\mathbf A) = \frac{1}{2\pi i} \oint_\gamma f(z)\, (z \mathbf I- \mathbf A)^{-1}\,\mathrm dz$$ where the ...

8

I needed this decomposition to answer another question, so I broke down and implemented it myself. The code is more or less a straightforward translation of the pseudocode in Golub/Van Loan: LDLT[mat_?SymmetricMatrixQ] := Module[{n = Length[mat], mt = mat, v, w}, Do[ If[j > 1, w = mt[[j, ;; j - 1]]; v = ...

7

$MachinePrecision is different from MachinePrecision. The former calls for an arbitrary precision calcluation, done at the same precision as the machine-precision one. The main reason one would want to use this is to enable precision tracking, which is absent for a true machine-precision calculation using MachinePrecision. And, there is your answer. ... 2 I write this answer with the caveat that I don't have Mathematica version 9 or later (the versions which now have this very belated function built-in), but with said caveat being offset by knowing a thing or two about the function of a matrix. ;) I'd have to agree with george's take that the docs for MatrixFunction[] could probably have explained things a ... 3 Here is a slightly compacted reformulation of belisarius's answer: a = Take[mat, 10, 10]; b = Take[mat, 10, -15]; c = Take[mat, -15, 10]; d = Take[mat, -15, -15]; rr = ArrayRules[d - c.SparseArray[LinearSolve[a, b]]]; detr = Det[SparseArray[rr /. (pos_ /; VectorQ[pos, IntegerQ] -> expr_) :> (pos -> C @@ pos), ... 0 A wasteful method: mat = {{2, 4, -1, 5, -2}, {-4, -5, 3, -8, 1}, {2, -5, -4, 1, 8}, {-6, 0, 7, -3, 1}}; {m, n} = Dimensions[mat]; {lu, piv} = Most[LUDecomposition[mat]] // Quiet; lm = Drop[LowerTriangularize[lu, -1], None, -1] + IdentityMatrix[m]; um = LinearSolve[lm, mat]; lm.um == IdentityMatrix[m][[piv]].mat True 4 Using as basis the great resource for core numerical algorithms below, I managed to implement a compiled linear solve which doesn't call MainEvaluate (so quite fast). I needed a linear solve for an optimization where the objective function requires inverting matrices, I was hesitating to use C++, but I preferred to stay in Mathematica. Resources ... 16 Using the properties of Block matrices: $$\det\begin{pmatrix}\mathbf A&\mathbf B\\\mathbf C&\mathbf D\end{pmatrix}=\det(\mathbf A)\det\left(\mathbf D-\mathbf C\mathbf A^{-1}\mathbf B\right)$$ To visualize your matrix: mat1 = mat; {mat1[[;; 10, ;; 10]], mat1[[;; 10, 11 ;;]], mat1[[11 ;;, ;; 10]], mat1[[11 ;;, 11 ;;]]} = Range@4; (* cool :) *) ... 2 Try this: (Normalize /@ evecs // Transpose // Assuming[0 < \[Theta] < Pi, FullSimplify[TrigToExp@#]] &) /. {(1 + Abs[Cot[\[Theta]/2]]^2)^(-1/2) -> Abs[Sin[\[Theta]/2]]} (* {{-I Cos[\[Theta]/2], I Sin[\[Theta]/2]}, {Sin[\[Theta]/2], Cos[\[Theta]/2]}} *) Edit: if theta<0, (Normalize /@ evecs // Transpose // ... 0 On the basis of$A$being defined as posted: adef[g_List, h_] := Module[{n = Length[g]}, gs = ArrayFlatten[Table[g, {n}]]; hs = ArrayFlatten[ ReplacePart[ConstantArray[0, {n, n}], {i_, i_} :> h]]; hs - gs] To illustrate: adefex[g_List, h_] := Module[{n = Length[g]}, gs = ArrayFlatten[Table[g, {n}]]; hs = ArrayFlatten[ ... 5 Here's my relatively compact implementation of Glynn's formula, which incorporates the Gray code optimization: SetAttributes[GrayCode, Listable]; GrayCode[n_Integer] := BitXor[n, BitShiftRight[n]] permanent[mat_?MatrixQ] /; Equal @@ Dimensions[mat] := Module[{b = 2^(Length[mat] - 1)}, PadRight[{}, b, {1, -1}].(Times @@@ ... 11 The Pitsianis-Van Loan algorithm turns out to be surprisingly easy to implement in Mathematica: nearestKroneckerProduct[mat_?MatrixQ, dim1_?VectorQ, dim2_?VectorQ] /; TrueQ[Dimensions[mat] == dim1 dim2] := Module[{bv, cv, sig}, {bv, sig, cv} = SingularValueDecomposition[Flatten[ Map[Composition[Flatten, Transpose], ... 0 Some people (see The ubiquitous Kronecker product by Van Loan) have worked on finding two matrices$A, B$of specified size whose tensor product$A \otimes B$is closest (in a norm) to a given (larger) matrix$C$. That is, find$A, B$which minimize$||C-A \otimes B||\$. The algorithm is based on the SVD. There is a matlab implementation somewhere. It would ...

9

There is another option, using the relatively new tensor capabilities of Mathematica. This is pretty much copied from another answer by jose, but I don't need any assumptions here: TensorExpand[KroneckerProduct[X, X] + KroneckerProduct[-X, X]] (* ==> 0 *) TensorExpand[KroneckerProduct[2 X, 3 Y]] (* ==> 6 KroneckerProduct[X, Y] *) There is a ...

3

How about: av = Array[Subscript[a, ##] &, {2}]; bv = Array[Subscript[b, ##] &, {2}]; KroneckerProduct[av, bv] + KroneckerProduct[-av, bv] {{0, 0}, {0, 0}}

1

Second update -- I should state in simplest terms the issue the OP is facing. The set-up. The equations for a given a are eqs = Table[Power[Sqrt[λ2/(μ2 c2)], n] p[0, n + 1] == β[n] p[0, 0] + Sum[α[n, k] p[0, k], {k, 0, a}] /. NumVal, {n, 0, a - 1}]; All the variables involved in eqs are given by vars = Table[p[0, n], {n, 0, a}]; (* starts ...

6

As @Guess who it is. states in the comments, an overdetermined linear problem can be solved using Mathematica's LeastSquares[] functionality. To input your above system of equations: a = {{1, -2, 2, 2, 0, -2}, {1, -2, 2, -2, 0, 2}, {1, -2, -2, -2, -4, -2}, {1, 2, 2, -2, -4, -2}, {1, 2, -2, -2, 0, 2}, {1, 2, -2, 2, 0, -2}, {1, -2, -2, 2, 4, ...

5

We will go by solving one equation at a time and generating the corresponding replacement rules. Beware of possible numerical instabilities. The following is the equivalence between the code in your edited example and my code on the previous incarnation of this answer. I think this is enough for you to use it. Please note that the only claim on the ...

0

(This answer is incomplete, but I don't know if I'll get time to finish it soon.) LUDecomposition nominally requires its input matrix to be square, as per the documentation. As far as I can see, the output in this case is correct except in the last row: {lu, a, b} = LUDecomposition[{{2, 4, -1, 5, -2}, {-4, -5, 3, -8, 1}, {2, -5, -4, 1, ...

4

If you call Orthogonalize at the end, you're orthogonalizing the eigenvectors in a different order (i.e. after sorting on eigenvalue, rather than before). Orthogonalizing the same list in a different order usually gives a different output. Orthogonalize[{{1., 2}, {1, 3}}] (* {{0.447214, 0.894427}, {-0.894427, 0.447214}} *) Orthogonalize[{{1, 3}, {1., 2}}] ...

1

Except in special cases, the problem seems overdetermined. Consider for simplicity m = {{m11, m12, m13}, {m21, m22, m23}, {m31, m32, m33}}; n = {{n11, n12, n13}, {n21, n22, n23}, {n31, n32, n33}}; t = {{t11, t12, t13}, {t21, t22, t23}, {t31, t32, t33}}; Then, the first equation becomes t.m - n.t == 0 unless t is singular. The solution is Solve[t.m - ...

Top 50 recent answers are included