2

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]
      ]
qahtah
  • 1,397
  • 6
  • 14

1 Answers1

2

[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}} *)`
Daniel Lichtblau
  • 58,970
  • 2
  • 101
  • 199
  • Thanks Daniel, actually the function is NsymDisp[5,w]. Reduce and NSolve don't give the same results. – qahtah Nov 13 '17 at 18:07
  • 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. – Daniel Lichtblau Nov 13 '17 at 19:50
  • 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. – qahtah Nov 15 '17 at 03:52