2
$\begingroup$

I am trying to find multiple complex and real roots of the determinant of symbolic matrix. Someone advised me to use Muller's method. I am wondering if there is a package already written to do this method?

The problem is: The data is given and the function to be solved for is NsymDisp[w_, k_], we solve NsymDisp for k if some values of w is given. 

Cnorm = 129.5;
hnorm = 5*10^-6;
data = {a -> 50*10^-6/hnorm,
 h -> 5*10^-6/hnorm,
  \[Mu] -> (84.7)/Cnorm,
  \[Lambda] -> (129.5)/Cnorm,
  Subscript[C, 11] -> (268)/Cnorm,
  Subscript[C, 13] -> (95)/Cnorm,
   Subscript[C, 33] -> (186)/Cnorm,
   Subscript[C, 55] -> (49)/Cnorm,
   Subscript[C, 44] -> (37)/Cnorm,
   Subscript[C, 66] -> (95)/Cnorm,
   Subscript[\[Rho], c] -> 1};
  NsymDisp[w_, k_] := 
  Module[{\[Omega] = w, Cp, Cs, p, 
   q, \[Alpha]11, \[Beta]11, \[Gamma]11, \[Alpha]12, \[Beta]12, \
 \[Gamma]12, \[Alpha]21, \[Beta]21, \[Gamma]21, \[Alpha]22, \[Beta]22, \
 \[Gamma]22},
  Cp = Sqrt[(\[Lambda] + 2 \[Mu])/Subscript[\[Rho], c]] /. data;
   Cs = Sqrt[\[Mu]/Subscript[\[Rho], c]] /. data;
   kp = \[Omega]/Cp;
   ks = \[Omega]/Cs;
   p = Sqrt[kp^2 - k^2];
   q = Sqrt[ks^2 - k^2];   
   \[Alpha]11 = 2 \[Mu] k p;
   \[Beta]11 = (Subscript[\[Rho], 
   c] \[Omega]^2 - (Subscript[C, 11] + 
     Subscript[C, 
      13] (\[Lambda] - Subscript[C, 13])/Subscript[C, 33]) k^2 - 
      Subscript[C, 13] (\[Lambda] + 2 \[Mu]) p^2/Subscript[C, 33]) k;
    \[Gamma]11 = (Subscript[C, 13]^2 (Subscript[C, 55] - 2 \[Mu]) k^2 - 
      Subscript[C, 
       33] (Subscript[C, 55] - 2 \[Mu]) (Subscript[C, 11] k^2 - 
         Subscript[\[Rho], c] \[Omega]^2) - 
      Subscript[C, 13] Subscript[C, 
       55] (2 \[Mu] k^2 - 
         Subscript[\[Rho], c] \[Omega]^2)) k p/(2 Subscript[C, 33]
        Subscript[C, 55]);
    NAs11 = ( \[Alpha]11 + h^2 \[Gamma]11)*Sin[p a] + 
     h \[Beta]11 Cos[p a];
    \[Alpha]12 = \[Mu] (q^2 - k^2);
    \[Beta]12 = (Subscript[\[Rho], 
        c] \[Omega]^2 - (Subscript[C, 11] - 
          Subscript[C, 
           13] (Subscript[C, 13] + 2 \[Mu])/Subscript[C, 33]) k^2) q;
     \[Gamma]12 = (Subscript[C, 
        13]^2 ((\[Mu] - Subscript[C, 55]) k^2 - \[Mu] q^2) k^2 + 
       Subscript[C, 
        33] ((Subscript[C, 55] - \[Mu]) k^2 + \[Mu] q^2) (Subscript[C, 
           11] k^2 - Subscript[\[Rho], c] \[Omega]^2) - 
        Subscript[C, 13] Subscript[C, 
        55] (Subscript[\[Rho], 
           c] \[Omega]^2 + \[Mu] (q^2 - k^2)) k^2)/(2 Subscript[C, 33]
         Subscript[C, 55]);
      NAs12 = (\[Alpha]12 + h^2 \[Gamma]12)*Sin[q a] + 
      h \[Beta]12 Cos[q a];          
     \[Alpha]21 = -\[Lambda] k^2 - (\[Lambda] + 2 \[Mu]) p^2;
     \[Beta]21 = (Subscript[\[Rho], c] \[Omega]^2 - 2 \[Mu] k^2) p;
       \[Gamma]21 = ((\[Lambda] + 
         2 \[Mu]) (Subscript[\[Rho], c] \[Omega]^2 + 
            Subscript[C, 13]
             k^2) p^2 - ((Subscript[C, 13] + Subscript[C, 
             33] - \[Lambda]) Subscript[\[Rho], 
            c] \[Omega]^2 - (Subscript[C, 11] Subscript[C, 33] + 
              Subscript[C, 
              13] (\[Lambda] - Subscript[C, 
                13])) k^2) k^2)/(2 Subscript[C, 33]);
      NAs21 = (\[Alpha]21 + h^2 \[Gamma]21) Cos[p a] + 
      h \[Beta]21 Sin[p a];
    \[Alpha]22 = 2 \[Mu] k q;
    \[Beta]22 = (\[Mu] (k^2 - q^2) - 
         Subscript[\[Rho], c] \[Omega]^2) k;
      \[Gamma]22 = ((Subscript[C, 11] Subscript[C, 33] - 
           Subscript[C, 
           13] (Subscript[C, 13] + 2 \[Mu])) k^2 - (Subscript[C, 13] + 
            Subscript[C, 33] + 2 \[Mu]) Subscript[\[Rho], 
        c] \[Omega]^2) k q/(2 Subscript[C, 33]);
   NAs22 = (\[Alpha]22 + h^2 \[Gamma]22 ) Cos[q a] + 
       h \[Beta]22 Sin[q a];             
     symMat = {{NAs11, NAs12}, {NAs21, NAs22}} /. data;
    Det[symMat // Chop]
      ]
Improve this question
$\endgroup$
7
  • 2
    $\begingroup$ may be this muller-method-in-mathematica $\endgroup$
    – Nasser
    Nov 12, 2017 at 4:30
  • 1
    $\begingroup$ If the "symbolic matrix" is not too complicated, can you include it here? $\endgroup$
    – J. M.'s eventual burnout ♦
    Nov 12, 2017 at 4:34
  • $\begingroup$ If you are satisfied to find the roots by any means, there are several routines described on this site. For instance, if the determinant is real, then question [16439] (mathematica.stackexchange.com/q/16439/1063) and its answers provide multiple approaches. $\endgroup$
    – bbgodfrey
    Nov 12, 2017 at 14:26
  • $\begingroup$ Absent a specific example, there is not really much of a question here. No clear indication of why Roots or NRoots would not suffice, for example. $\endgroup$
    – Daniel Lichtblau
    Nov 12, 2017 at 15:30
  • 1
    $\begingroup$ Apparently i had to be validated as "Not a robot". Now I can add comments. my hat's off to MSE for such an informative indication (read: none whatsoever) of what was the problem. $\endgroup$
    – Daniel Lichtblau
    Nov 12, 2017 at 17:23

1 Answer 1

Reset to default
2
$\begingroup$

[Not a full answer but I m having browser problems and cannot add a comment.]

If a range is specified then NSolve can find roots in that range.

In[22]:= Chop[NSolve[NsymDisp[w, 5] == 0 && Abs[w] < 6, w]]

(* Out[22]= {{w -> -5.83531}, {w -> -5.65654}, {w -> -5.31554}, {w -> \
-5.1549}, {w -> -5.00189}, {w -> -4.85735}, {w -> -4.72215}, {w -> \
-4.59719}, {w -> -4.48334}, {w -> -4.38144}, {w -> -4.29225}, {w -> \
-4.21635}, {w -> -4.15414}, {w -> -4.10577}, {w -> -4.07127}, {w -> \
-4.05058}, {w -> -4.04368}, {w -> 4.05058}, {w -> 4.07127}, {w -> 
4.10577}, {w -> 4.15414}, {w -> 4.21635}, {w -> 4.29225}, {w -> 
4.38144}, {w -> 4.59719}, {w -> 4.72215}, {w -> 5.00189}, {w -> 
5.1549}, {w -> 5.31554}, {w -> 5.65654}, {w -> 5.83531}} *)`
Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks Daniel, actually the function is NsymDisp[5,w]. Reduce and NSolve don't give the same results. $\endgroup$
    – qahtah
    Nov 13, 2017 at 18:07
  • $\begingroup$ Are you sure they are different? NSolve[NsymDisp[5, w] == 0 && Abs[w] < 3, w] and Reduce[NsymDisp[5, w] == 0 && Abs[w] < 3, w] look likew the same results to my eye. This is using Mathematica 11.2, in case that matters. $\endgroup$
    – Daniel Lichtblau
    Nov 13, 2017 at 19:50
  • $\begingroup$ Thanks Daniel, I am not talking about a single point, I varied the value of w from 0 to 10 with an increment of 0.1. $\endgroup$
    – qahtah
    Nov 15, 2017 at 3:52

Your Answer

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

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