0
$\begingroup$

When I try to solve my question according to @DanielLichtblau's comments, I encountered another issue, with which I have been struggling for a whole night.

Taking a matrix for example,

mat[a_, b_, c_, k_] = {{Cos[a*k], f^3 + 2, -I*f + k/4}, {0, b*Sin[a*k], f}, {1/f^2, (k^2 - 1)/c, k*Cos[a*k]}};

where a, b, c and k are all real numbers, while f is a complex variable

Defining its determinant with a given parameter c=1

det[a_, b_, k_] := Det[mat[a, b, 1, k]]

Because I need the determinant to be zero, taking its numerator should be enough. Also, I assign values to the remaining parameters (a, b, k).

fexpr = Numerator[det[1, 1, 0.1]]

In view of the variable f in the matrix being complex-valued, expanding the numerator and making explicit replacement of its real and imaginary parts

polys = ComplexExpand[ReIm[fexpr], f] /. Re[f] -> ref /. Im[f] -> imf

Solving for the real and imaginary parts of f that make the numerator of the determinant vanishing

NSolve[polys, {imf, ref}, Reals]

Here the constraint of Reals was used since both imf and ref are real-valued, we only need the real-valued solutions. But, it reports an error which appears to relate to the Reals constraint.

NSolve::nddc: The system (-0.00249584+0. I)-0.0998334 imf-0.00988384 imf^2+1. imf^4+2. ref-2.95516 imf^2 ref+0.00988384 ref^2-6. imf^2 ref^2+0.985054 ref^3+1. ref^4==0&&(0. +0. I)+<<7>>+2.95516 imf ref^2+4. imf ref^3==0 contains a nonreal constant -0.00249584+0. I. With the domain [DoubleStruckCapitalR] specified, all constants should be real.

Strangely, if I try a slightly modified polys1

polys1 = {-0.1 - 0.1*imf - 0.01*imf^2 + 1.1*imf^4 + 2.01*ref - 
3.9*imf^2 ref + 0.01*ref^2 - 6.1*imf^2 ref^2 + 0.98*ref^3 + 
1.01*ref^4, 0.01 + 2.01*imf - 0.985*imf^3 + 0.02*imf ref - 4.01*imf^3 ref + 2.955*imf ref^2 + 4*imf ref^3};

It turns to work well with NSolve

NSolve[polys1, {imf, ref}, Reals]

(*{{imf -> 0.00113843, ref -> -1.68558}, {imf -> 1.05537, ref -> 0.278329}, {imf -> -1.03992, ref -> 0.310137}, {imf -> -0.00495377, ref -> 0.0494334}}*)

Can anyone please help with this. Thank you very much!

Update (Dec22)

I am trying to combine @Bob's answer into a loop but the modified code always gave many errors. I did slight modifications, which can give correct results without looping. Could you give any suggestion?

det[a_, b_, k_] = With[{a = a, b = b, k = k}, Det[mat[a, b, 1, k]]];

CurvePoints = {};

fr[k_?NumericQ, a_?NumericQ, b_?NumericQ] := Module[{sol, solrAbs, fexpr}, 
fexpr = Numerator[ComplexExpand[det[a, b, k], {f}]] /. f -> fre + I*fim;
sol = {f -> fre + I*fim} /. NSolve[Thread[(ReIm[fexpr] // ComplexExpand) == 0], {fre, fim}, Reals, WorkingPrecision -> 20];
solrAbs = Abs[Re@sol[[All, 1, 2]]];
Table[If[solrAbs[[n]] < 0.35,
AppendTo[CurvePoints, {k, a, b} -> sol[[n, 1, 2]]]; 
fr[k, a, b] = {k, a, b, sol[[n, 1, 2]]}, Sequence[]], {n, 1, Length[sol]}]]

This is the loop with a temporary monitor

pts = Monitor[Table[fr[k, a, b], {a, 0, 1, 0.5}, {b, 0, 1, 0.5}, {k, 0, 0.1, 0.1}], {k, a, b, s}];
$\endgroup$
4
  • $\begingroup$ Are you looking for an analytical solution? If not what are the values of parameter a,b,k? $\endgroup$ Dec 21, 2021 at 10:05
  • $\begingroup$ @UlrichNeumann I am looking for numerical solutions. Please see my update. $\endgroup$
    – Enter
    Dec 21, 2021 at 10:08
  • $\begingroup$ Why do you not simply write: Solve[fexpr == 0, f] ? $\endgroup$ Dec 21, 2021 at 10:28
  • $\begingroup$ @DanielHuber ah! The reason is that I need to handle the real and imaginary parts of f separately in my real problem. Please see my previous question and discussion therein. $\endgroup$
    – Enter
    Dec 21, 2021 at 11:22

2 Answers 2

2
$\begingroup$
mat[a_, b_, c_, 
   k_] = {{Cos[a*k], f^3 + 2, -I*f + k/4}, {0, b*Sin[a*k], 
    f}, {1/f^2, (k^2 - 1)/c, k*Cos[a*k]}};

det[a_, b_, k_] = 
  ComplexExpand[Det[mat[a, b, 1, k]], {f}] // 
   FullSimplify;

fexpr = Numerator[det[1, 1, 1/10]] /. 
  f -> fre + I*fim;

sol = {f -> fre + I*fim} /. 
   NSolve[
    Thread[(ReIm[fexpr] // ComplexExpand // 
     FullSimplify) == 0], {fre,
      fim}, Reals, WorkingPrecision -> 20] // N

(* {{f -> 0.33036 + 1.04213 I}, 
    {f -> -1.68481 - 0.0191733 I}, 
    {f -> 0.368153 - 1.02289 I}, 
    {f -> 0.00124481 - 0.0000621358 I}} *)

det[1, 1, 1/10] /. sol

(* {-9.83382*10^-17 - 1.82962*10^-16 I, 
     4.4498*10^-19 + 1.95483*10^-17 I,
     5.98657*10^-17 + 7.23935*10^-17 I, 
     1.73762*10^-15 - 1.7362*10^-14 I} *)

EDIT: For your revised question. The imaginary terms that were occurring in the NSolve were artifacts from the use of machine precision.

mat[a_, b_, c_, 
   k_] = {{Cos[a*k], f^3 + 2, -I*f + k/4}, {0, b*Sin[a*k], 
    f}, {1/f^2, (k^2 - 1)/c, k*Cos[a*k]}};

det[a_, b_, k_] = With[{a = a, b = b, k = k}, Det[mat[a, b, 1, k]]];

CurvePoints = {};

fr[k_?NumericQ, a_?NumericQ, b_?NumericQ] := 
 Module[{sol, solrAbs, fexpr}, 
  fexpr = Numerator[FullSimplify[
      ComplexExpand[det @@ SetPrecision[{a, b, k}, 35], {f}]]] /. 
    f -> fre + I*fim;
  sol = {f -> fre + I*fim} /. NSolve[
     Thread[(ReIm[fexpr] // ComplexExpand // FullSimplify) == 0], 
     {fre, fim}, Reals, WorkingPrecision -> 20];
  solrAbs = Abs[Re@sol[[All, 1, 2]]];
  Table[If[solrAbs[[n]] < 7/20, 
    AppendTo[CurvePoints, {k, a, b} -> sol[[n, 1, 2]]];
    fr[k, a, b] = {k, a, b, sol[[n, 1, 2]]}, Sequence[]],
   {n, 1, Length[sol]}]]

Off[NSolve::precw]

pts = Monitor[
   Table[fr[k, a, b], {a, 0, 1, 1/2}, {b, 0, 1, 1/2}, 
     {k, 0, 1/10, 1/10}], {k, a, b, s}];

The format of the CurvePoints is unusual

CurvePoints[[1]]

(* {0, 0, 0} -> 0.3478103847799310287 + 1.0288522541366924034  I *)
$\endgroup$
6
  • $\begingroup$ Thank you for the answer! I am trying to modify your code into a loop but the modified code always gave errors. I did some slight modifications only, which can give the same result as yours without looping. Could you give some suggestion? $\endgroup$
    – Enter
    Dec 22, 2021 at 6:19
  • $\begingroup$ impressive! It now works at least. But checking pts, it is found some elements are something like If[Re[{f -> I fim + fre}[[All, 1, 2]]]...] instead of numeric. Also each of the unevaluated element seems to correspond to a warning of Part specification {f->I\ fim+fre}[[All,1,2]] is longer than depth of object. Why pts include these non-numeric elements? Does it mean NSolve fail to find a solution for the corresponding parameter values. If this is the case, is there any way to remove them? $\endgroup$
    – Enter
    Dec 22, 2021 at 9:12
  • $\begingroup$ I have checked with the parameter values of the non-numeric elements and found that the sols exist actually with those parameter values ... $\endgroup$
    – Enter
    Dec 22, 2021 at 9:29
  • $\begingroup$ I do not see any error messages nor any If statements in pts. There are Null from the Sequence[]. Change Sequence[] to Nothing to remove the Null $\endgroup$
    – Bob Hanlon
    Dec 22, 2021 at 16:48
  • $\begingroup$ Thread[(ReIm[fexpr] // FullSimplify) == 0] by itself does not ensure that all of the complex terms have been cancelled out. The ComplexExpand helps to ensure that the necessary transformations are performed to eliminate any complex representations. If any complex representations remain, the NSolve will fail since its domain is restricted to Reals $\endgroup$
    – Bob Hanlon
    Dec 24, 2021 at 4:24
1
$\begingroup$

Try NMinimize

mat[a_, b_, c_,k_] := {{Cos[a*k], f^3 + 2, -I*f + k/4}, {0, b*Sin[a*k],f}, {1/f^2, (k^2 - 1)/c, k*Cos[a*k]}}
/. f -> fre + I fim;
det=Det[mat[1, 1, 1, 1/10]] // Simplify
eqn = Thread[ComplexExpand[ReIm[det]] == 0]

NMinimize evaluates a minimum solution which fullfills the equations

sol=NMinimize[{fre^2 + fim^2, eqn}, {fre, fim}]
(*{1.19517, {fre -> 0.33036, fim -> 1.04213}}*)

eqn /. Equal -> Subtract /. sol[[2]]
(*{2.14508*10^-9, -4.62391*10^-8}*)

It seems to be a unique solution

Show[{ContourPlot[Evaluate [eqn], {fre, -.020, 1   }, {fim, 0, 1.5 }]  ,
Graphics[Point[{fre,fim}/. sol[[2]]]]}] 

enter image description here

Hope it helps!

addendum

NSolveworks too inside a given parameter range:

NSolve[{eqn, 0 < fre < 2, 0 < fim < 2} // Flatten, {fre, fim}]
(*{{fre -> 0.33036, fim -> 1.04213}}*)
$\endgroup$
4
  • $\begingroup$ It's only necessary if the solution isn't unique. $\endgroup$ Dec 21, 2021 at 13:37
  • $\begingroup$ in your case also NMinimize[{1, eqn}, {fre, fim}] or NMinimize[{ #.#&[eqn/.Equal->Subtract]}, {fre, fim}] give the same result! $\endgroup$ Dec 21, 2021 at 13:44
  • $\begingroup$ What means "What is the rationale?"? $\endgroup$ Dec 21, 2021 at 13:46
  • $\begingroup$ now I understand your method, but I think Bob's method is more suitable for my problem. Thank you again! $\endgroup$
    – Enter
    Dec 21, 2021 at 14:18

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.