0
$\begingroup$

I'm trying to use Mathematica to find the numerical solution to an equation Det[M[x]] == 0, where M[x] is a matrix function of x, defined below:

r = 10; ω = E^(π I/r); 
M[x_] = Table[D[Cos[ω^j x], {x, i}], {i, 0, r - 1}, {j, 0, r - 1}]; 
G[x_] := Det[N[M[x]]];
M[x] // MatrixForm

I expect to see the solution to G[x] == 0 for large r up to a few hundreds, but when I call the function FindRoot, it is already very slow even for r = 10 (doesn't finish in a few seconds):

FindRoot[G[x] == 0, {x, (r + 1)/4 π}]

But the plot of G[x] is extremely fast:

Plot[Norm[G[x]], {x, (r + 0.6)/4 π, (r + 0.7)/4 π}]

And actually, I tried manually using bisection method and plotting over and over again (each time decreasing the interval by a factor of 2), which give me a very accurate solution in just a minute or so, much faster than FindRoot.

So, why is FindRoot so slow in this case?

$\endgroup$

2 Answers 2

2
$\begingroup$

I think the problem here is that you're using Det to detect if the matrix is singular. Determinants are very messy and erratic functions and are generally not a very good measure for how close a matrix is to being singular. That's why FindRoot struggles. Instead, you can search for a matrix for which the least singular value is smallest:

ClearAll[M, G1, G2]
r = 10;
ω = E^(\[Pi] I/r);
M[x_] = Table[D[Cos[ω^j x], {x, i}], {i, 0, r - 1}, {j, 0, r - 1}];
G1[x_?NumericQ] := Det[M[x]];
G2[x_?NumericQ] := First @ SingularValueList[M[x], -1];

sol = NMinimize[G2[x], x, WorkingPrecision -> 50]

{0, {x -> 0.46507640237572384522369322963792665762298255262859}}

Check solution:

N @ G1[x /. Last[sol]]

2.4946*10^-40 + 3.6561*10^-90 I

Edit I just noticed that I left the N function inside of G2, which is probably not a good idea. Also, you can use FindMinimum instead of NMinimize to get your answer a little faster. It's not unique, obviously.

$\endgroup$
3
  • $\begingroup$ Thanks for your answer. I'm also interested in the null vector of M[x] at the point where it is singular. I tried NullSpace[M[x0]] (where x0 is the solution to the equation G[x]==0), but it gives nothing, probably because numerical solution is only approximate. So is there an easy way to find this null vector by slightly modifying your code? $\endgroup$
    – Lagrenge
    Commented Oct 5, 2019 at 16:39
  • $\begingroup$ @Lagrenge I'd have to think about that for a bit, but you probably want to read up on SingularValueDecomposition. That sounds more or less like what you would need, I think, though I'd have to piece together the details for myself as well. I also believe that NullSpace has a Tolerance option. You may want to try that as well. $\endgroup$ Commented Oct 5, 2019 at 19:32
  • $\begingroup$ @Lagrenge This answer might be of use to you: math.stackexchange.com/questions/1771013/… $\endgroup$ Commented Oct 5, 2019 at 19:37
1
$\begingroup$
Clear[G];
G[x_?NumericQ] := Det[N[M[x]]]

Then FindRoot is not too slow:

FindRoot[G[x] == 0, {x, (r+1)/4 Pi}] //AbsoluteTiming

FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances.

{0.017618, {x -> 8.44175 + 5.15677*10^-13 I}}

$\endgroup$
3
  • $\begingroup$ It seems like the FindRoot's result is a very poor solution though, using those numbers G[x] returns -196217.-147916 I. WorkingPrecision doesn't improve the accuracy either. $\endgroup$
    – eyorble
    Commented Oct 5, 2019 at 5:02
  • $\begingroup$ It seems accuracy improves significantly if the call to N is removed from G[x]'s definition. I would think it's redundant once ?NumericQ is applied anyways. $\endgroup$
    – eyorble
    Commented Oct 5, 2019 at 5:22
  • $\begingroup$ @eyorble It is not redundant if Det and N are far from commuting. I have not checked for this example but it would not surprise me, especially if either of Det or FindRoot is using Expand. $\endgroup$ Commented Oct 5, 2019 at 15:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.