1
$\begingroup$

I want to solve the following problem: Given a complex matrix $m=\textbf{m}(J, ω)$, with $J, ω \in \mathbb R$, find $λ$ and $ω$ such that

$$\det[\textbf{m} (J, ω)] = 0$$

I tried to solve it with FindRoot. I've found a solution with specific initial values $(J_c,\omega_c)$, but it works only with MachinePrecision->13 and not grater MachinePrecision. Moreover, when I test the solution, it gives me back a zero with a low accuracy, about $10^-3$.

My question is the following: is it possible that FindRoot finds points which are not actual zeros of the function but point for which the function is very small?

How can I distinguish a zero for a "almost zero" numerically?

d = Rationalize[5/10^6]; \[Tau] = Rationalize[0.1]; L = Rationalize[15/10^3]; h = Rationalize[50/10^6]; \[Xi] = Rationalize[0.001]; \[Rho] = Rationalize[1400]; 
  \[CapitalGamma] = Rationalize[0.001]; c = Rationalize[0.05]; Ee = Rationalize[10^8]; w = Rationalize[10^(-3)]; 
A = Rationalize[w*h]; \[CapitalIota] = Rationalize[w*(h^3/12)]; B = Rationalize[Ee*\[CapitalIota]]; \[Alpha] = Rationalize[(d/2)*(2*d - h - Exp[-h/d]*(2*d + h))]; 
  \[Lambda] = Rationalize[c*\[CapitalGamma]*\[Alpha]*(12/h^3)]; 

m2[\[Beta]j_, \[CapitalLambda]_, \[Sigma]_] = ( \[Beta]j^3 - \
\[CapitalLambda] \[Beta]j^2)
m3[\[Beta]j_, \[CapitalLambda]_, \[Sigma]_] = ( \[Beta]j^2 - \
\[CapitalLambda] \[Beta]j) Exp[\[Beta]j L]; 
m4[\[Beta]j_, \[CapitalLambda]_, \[Sigma]_] = ( \[Beta]j - \
\[CapitalLambda]) Exp[\[Beta]j L]; 

(*To find \[Beta]j (j=1,2,3,4) we solve a polynomial equation.*)
eigen[\[CapitalLambda]_, \[Sigma]_] = 
  FullSimplify[
   B ( x^4 - \[CapitalLambda] x^3) + \[Sigma] (\[Xi] + \[Rho] A \
\[Sigma])];

root = Simplify[Solve[eigen[\[CapitalLambda], \[Sigma]] == 0, x]];
xx[\[CapitalLambda]_, \[Sigma]_] = x /. root;

\[CapitalLambda]lin[J_, \[Sigma]_] = 
 FullSimplify[\[Lambda] J /(1/\[Tau] + \[Sigma])]

m[\[CapitalLambda]_, \[Sigma]_] = {{m2[
     xx[\[CapitalLambda], \[Sigma]][[1]], \[CapitalLambda], \[Sigma]],
     m2[xx[\[CapitalLambda], \[Sigma]][[
      2]], \[CapitalLambda], \[Sigma]], 
    m2[xx[\[CapitalLambda], \[Sigma]][[
      3]], \[CapitalLambda], \[Sigma]], 
    m2[xx[\[CapitalLambda], \[Sigma]][[
      4]], \[CapitalLambda], \[Sigma]]}, {1, 1, 1, 
    1}, {m3[xx[\[CapitalLambda], \[Sigma]][[
      1]], \[CapitalLambda], \[Sigma]], 
    m3[xx[\[CapitalLambda], \[Sigma]][[
      2]], \[CapitalLambda], \[Sigma]], 
    m3[xx[\[CapitalLambda], \[Sigma]][[
      3]], \[CapitalLambda], \[Sigma]], 
    m3[xx[\[CapitalLambda], \[Sigma]][[
      4]], \[CapitalLambda], \[Sigma]]}, {m4[
     xx[\[CapitalLambda], \[Sigma]][[1]], \[CapitalLambda], \[Sigma]],
     m4[xx[\[CapitalLambda], \[Sigma]][[
      2]], \[CapitalLambda], \[Sigma]], 
    m4[xx[\[CapitalLambda], \[Sigma]][[
      3]], \[CapitalLambda], \[Sigma]], 
    m4[xx[\[CapitalLambda], \[Sigma]][[
      4]], \[CapitalLambda], \[Sigma]]}};



detexp[J_, \[Sigma]_] = 
  Det[m[\[CapitalLambda]lin[J, \[Sigma]], \[Sigma]]]/10^12;

det[J_, \[Omega]_] = detexp[J, I \[Omega]];

Jc = 4.5*10^3 // Rationalize;
\[Omega]c = 15;

firstmode = FindRoot[{Re[det[J, \[Omega]]] == 0, Im[det[J, \[Omega]]] == 0}, {J, 
  Jc}, {\[Omega], \[Omega]c}, WorkingPrecision -> 13]

(*out {J -> 4397.165861009, \[Omega] -> 14.81874398305}*)

{J1, \[Omega]1} = {J, \[Omega]} /. firstmode;

det[J1, \[Omega]1]

(*out 0.*10^-3 + 0.*10^-3 I*)

PS. I asked a similar question How can I improve the accuracy of FindRoot for a complex equation? and tried to implement the old suggestions, but the question here is slightly different.

Thank you!

Improve this question
$\endgroup$
3
  • $\begingroup$ First: You cannot have _ like \[Omega]_c in variable names in MMA. Second: Rationalize does noting to values that are already rational and only makes people question that and other things in your code. Third {Re[det[J, \[Omega]]]==0,Im[det[J, \[Omega]]]==0} seems like det[J,\[Omega]]]==0 Fourth: If you want an exact zero then I'd try eliminating all decimal points and possibly simplifying the problem a little and finally ask Reduce to find the exact solutions and see what happens. $\endgroup$
    – Bill
    Commented Jul 24 at 19:55
  • $\begingroup$ Hi, thank you for your answer. The first thing is now correct. The second thing: I see what you mean, but it should not affect the result. The third thing: the two equations are needed since we have two unknowns $(J,\omega)$, and hence two equations to solve. Since the determinant is complex, I took the real and the imaginary part of it. $\endgroup$
    – fritess
    Commented Jul 24 at 20:39
  • 1
    $\begingroup$ Why don't you like the answers you get with higher WorkingPrecision? They certainly get the determinant closer to zero. $\endgroup$
    – Bill Watts
    Commented Jul 24 at 20:55

1 Answer 1

Reset to default
2
$\begingroup$

I didn't get the exactly same result as you report, but I got something that needed fixing. It seems a WorkingPrecision needs to be much higher than the AccuracyGoal for this problem.

firstmode = 
  FindRoot[{Re[det[J, \[Omega]]] == 0, Im[det[J, \[Omega]]] == 0}, {J,
     Jc}, {\[Omega], \[Omega]c}, AccuracyGoal -> 16, 
   WorkingPrecision -> 64];
{J1, \[Omega]1} = {J, \[Omega]} /. firstmode;

N@det[J1, \[Omega]1]
det[N@J1, N@\[Omega]1]
(*
-2.89809*10^-18 - 8.48501*10^-17 I
-3.41099*10^-8 - 1.30445*10^-8 I
*)

The last line shows the loss of accuracy when the high-precision result is rounded to machine precision.

Improve this answer
$\endgroup$

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.