# Eigenvectors are divided by zero depending on evaluation, 6x6 matrix

I calculated the Eigenvectors of the $$6\times 6$$ matrix $$m$$ with parameters $$(a,b,c,d,e,f)\in \mathbb{R_{\ge 0}}$$. If I set e.g. $$b\rightarrow 0$$ after calculation, then 4 of the 6 Eigenvectors are affected by division by zero. If I set $$b\rightarrow 0$$ before I calculate the Eigenvectors then division by zero does not take place. To set $$b\rightarrow 0$$ is only one of many situations when the denominator becomes zero, therefore I am not looking for a special solution that works only for this case. How can I avoid division by zero if I want to replace the values $$a,b,c,d,e,f$$ after I calculated the Eigenvectors?

This is a follow up question from this post. There we have a 4x4 matrix.

Code:

m={{0,0,0,a c,0,-a b},{0,0,-d e,0,e f,0},{0,-d e,0,0,0,b d},{a c,0,0,0,-c f,0},{0,e f,0,-c f,0,0},{-a b,0,b d,0,0,0}};
v=FullSimplify[Eigenvectors[m],{a>=0,b>=0,c>=0,d>=0,e>=0,f>=0}];
v/.b->0 (*ComplexInfinity*)
Limit[v, b -> 0, Direction -> "FromAbove"] (*DirectedInfinity*)
v=Eigenvectors[m/.b->0] (*valid vectors*)


MMA 13.3

Comment to the answers, Nov. 13th, 2023

The answers use tricks to avoid division by zero. But we learn nothing from this and don't see which steps take place. Actually I thought the expression for the Eigenvectors can be reformulated so that any division by zero is avoided at all. Considering the answers it would be easier to set the parameters to zero before we calculate the Eigenvectors. Then we also learn nothing but it is faster and without tricks.

• Clear denominators and see if that helps. denoms = Map[Apply[PolynomialLCM, Denominator[Together@#]] &, v]; newv = Table[Together[denoms[[j]]*v[[j]]], {j, 6}]; newv /. b -> 0 Commented Nov 7, 2023 at 21:56
• I've noticed that there are two eigenvectors with eigenvalue zero. The eigenvectors have very simple forms, and so the limit $b\to0$ is very simple. For the other four eigenvectors, the eigenvalues come in pairs: $\pm\lambda_1$ and $\pm\lambda_2$, where $\lambda_1\neq\lambda_2$ even when $b=0$. This non-degeneracy in the $b=0$ limit means that you can just compute the eigenvectors by setting $b=0$ first, and there won't be any problem, i.e., the eigenvectors computed with $b\neq0$ will smoothly connect to the eigenvectors computed after setting $b=0$ (barring branch cut issues). Commented Nov 7, 2023 at 22:04
• ...or possible overall negative signs between corresponding eigenvectors in the $b=0$ limit and $b=0$ case, which you can fix by hand. The point being, why do you need to compute the $b\neq0$ case and then take the limit when you can just compute $b\neq0$ and $b=0$ cases separately? Commented Nov 7, 2023 at 22:06
• $b=0$ is only one special case. Other variables and many combinations of them can also make the denominator zero. I want a general solution for the Eigenvectors that works for any value. There are many special cases and I do not want to calculate them separately. Commented Nov 7, 2023 at 22:24
• @DanielLichtblau Works for the special case in the OP if $b\rightarrow 0$ but fails for example for $b\rightarrow 0,c\rightarrow 0$. The Eigenvectors become null vectors. There are many combinations when the denominator is zero. So there should be a general solution. Commented Nov 7, 2023 at 22:46

## 4 Answers

This is not an answer but too long for a comment. I can try to explain why you might not want to go the route you seem to be taking.

First observe that, independent of the parameters, 0 is an eigenvalue with multiplicity 2. So you get two eigenvectors just by taking the null space.

NullSpace[m]

(* Out[28]= {{0, b/e, 0, b/c, 0, 1}, {f/a, 0, f/d, 0, 1, 0}} *)


Of course you can clear denominators. That's as I mentioned in a comment. You responded, correctly, that this has problems of its own, if, say, parameters b and c are both zero, because then an eigenvector vanishes. Also the multiplicity of the zero eigenvalue changes in that case, jumping to 4. All troublesome.

So one might instead use Reduce to get a breakdown by cases based on the parameters. It could be done as below.

vec = Array[v, 6];
rr = Reduce[m . vec == 0, Variables[{vec, m}]];


But this is a nontrivial disjunction.

In[49]:= {Head@rr, Length@rr}

(* Out[49]= {Or, 206} *)


And there are two other conjugate pairs of eigenvalues to consider, each with complicated parametrized radicals.

My best guess as to what might work is to compute denominators of your eigenvectors, and use the vanishing of those to break into subcases. For each subcase it will still be a situation of rinse-and-repeat though; that's basically what happened with your simple example of setting not just b but also c to zero.

Bottom line: This process can get into considerable complexity quite quickly.

This is too long for a comment, but not quite a full solution, since it would require tweaking for the particular problem at hand. But I can at least illustrate a method of fixing this issue that works for 2-by-2 symmetric matrices. Note:

m = {{a, b}, {b, -a}};
Limit[Eigenvectors@m, b -> 0]
(* {{Indeterminate, 1}, {Indeterminate, 1}} *)


so we have the same problem as in OP. But if we normalize the vectors, then we fix part of the problem:

eigs = Normalize /@ Eigenvectors@m // ComplexExpand // Simplify;
Limit[eigs, b -> 0, Assumptions -> {a > 0}]
(* {{0, 1}, {Indeterminate, 0}} *)


The Indeterminate part comes from the fact that the two-sided limit doesn't exist, which is a consequence that the limit is actually different depending on whether b is negative or positive. However, the one-sided limit do exist, and they yield the correct limits for the eigenvectors:

Limit[eigs, b -> 0, Assumptions -> {a > 0}, Direction -> 1]
(* {{0, 1}, {-1, 0}} *)


and

Limit[eigs, b -> 0, Assumptions -> {a > 0}, Direction -> -1]
(* {{0, 1}, {1, 0}} *)


The issue with applying this method to the matrix in the OP is that it might take a huge amount of computing time to do the symbolic algebra required to normalize, simplify, and take the limit of the eigenvectors. Nonetheless, the idea should work.

• Is eternally busy by calculating. Therefore not applicable. Commented Nov 7, 2023 at 3:02
• @granularbastard I suspect that no automated method is going to work, then. I think that you have to normalize the vectors in order to render every element of the vector finite in the limit, and this is just going to be hard to do symbolically in general. However, there might be some by-hand simplifications one can do with your particular matrix. Given that this is the second question you've asked like this, I suspect you might have multiple matrices you want to do this with, so a general method won't work, but would you accept one that works for this special case only? Commented Nov 7, 2023 at 3:12
• Here we have the real problem. The first post (mathematica.stackexchange.com/q/292408/69288) was a simplified version, as I thought it is easier to work with a simpler example. However the solution for the simple version fails for the real problem. Therefore I had to open a new post. Commented Nov 7, 2023 at 17:47

This is just an application of answer posted here.

m = {{0, 0, 0, a  c, 0, -a  b}, {0, 0, -d  e, 0, e  f, 0}, {0, -d  e,
0, 0, 0, b  d}, {a  c, 0, 0, 0, -c  f, 0}, {0, e  f, 0, -c  f, 0,
0}, {-a  b, 0, b  d, 0, 0, 0}};
$Assumptions = a >= 0 && b >= 0 && c >= 0 && d >= 0 && e >= 0 && f >= 0; (*b->0 after calculation of Eigensystem*) v = FullSimplify[Eigenvectors[m]] Limit[(#/ Fold[If[Abs[ Limit[#2/#1, b -> 0, Direction -> "FromAbove"] /. Thread[{a, b, c, d, e, f} -> RandomReal[{1, 2}, 6]]] < Infinity, #1, #2] &, #]), b -> 0, Direction -> "FromAbove"] & /@ v  Basically what we are doing here is first use Fold to find out the maximum component as b->0 and then divide the whole eigenvector by that component to get an eigenvector with elements not approaching infinity as b->0. • If e.g. I try simultaneously$b\rightarrow0,c\rightarrow0\$ then I get messages like "Limit::ivar: 1.608919538582326 is not a valid variable.". So it seems this method is not very reliable. Commented Nov 13, 2023 at 0:39
• @granularbastard This is because some limits are not well-defined when both b and c approaches 0, and in such cases Mathematica could refuse to compute its limit and throw out the expression instead. To fix this simply change Thread[{a, b, c, d, e, f} -> RandomReal[{1, 2}, 6]] to Thread[{a, d, e, f} -> RandomReal[{1, 2}, 4]] and it works. Or if the order of limit don't matter, just change the way you take the limit to something like Limit[Limit[#2/#1, b -> 0, Direction -> "FromAbove"], c -> 0, Direction -> "FromAbove"]
– Wjx
Commented Nov 13, 2023 at 8:24

Eigenvectors can be expanded by multiplying with an arbitrary number. So why not just multiply each eigenvector with the product of the denominators of its components?

m = {{0, 0, 0, a c, 0, -a b}, {0, 0, -d e, 0, e f, 0}, {0, -d e, 0, 0,
0, b d}, {a c, 0, 0, 0, -c f, 0}, {0, e f, 0, -c f, 0, 0}, {-a b,
0, b d, 0, 0, 0}};
es = Eigensystem[m];
eigenvectors =  FullSimplify[Table[es[[2, i]] Times @@ Map[Denominator, es[[2, i, ;;]]], {i, 1, 6}]]


The results are indeed the eigenvectors:

Table[(m - es[[1, i]] IdentityMatrix[6]) . eigenvectors[[i]] // FullSimplify, {i, 1, 6}]


but the components don't contain any fractions. I have given more details on what is going on here.

• After eigenvectors /.b->0` all the eigenvectors corresponding to non-zero eigenvalues are zero vectors. Commented Nov 24, 2023 at 9:10