3

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

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

Since one can take the real and imaginary parts of the determinant, we end up with two equations.

I tried to solve it with FindRoot. Mathematica seems to find the right solutions for λ and ω, but when I try to plug them in the determinant, the determinant gives back a value of order $10^{-10}$.

I would like the solution to be more accurate so that the determinant would be of order $10^{-12}$ or less, is it possible?

The code is the following:

(* Parameters *)
Clear[sp]
ll = 58.3;
xin = 0.0034;
b = 1730;
a = 0.185;
k = 228;
e = 0.1;
t = 0.0038;
l0 = ( b t / xin)^(1/4);
ln[sp_] = sp*l0;
oursp = ll/l0

(* Compute the elements of the 4x4 matrix m.  This procedure is very specific to the problem that I’m working with.  *)


(bc1, bc2, bc4 will be of m)
bc1[xj_, c_] = (b xj^3 - c a^2 xj);
bc2[xj_] = Exp[xj ln[sp]] xj;
bc4[xj_, c_] = (Exp[xj ln[sp]]) (b xj^2 - c a^2 );


(To find xj (j=1,2,3,4) we solve a polynomial equation.)
eigen[c_, s_] = b x^4 - a^2 c x^2 + xin  s;
root = NSolve[eigen[c, s] == 0, x]; 
xx[c_, s_] = x /. root ;


(* Compute the 4x4 matrix and the determinant. *)
m[c_, s_, 
   sp_] = {{bc1[xx[c, s][[1]], c], bc1[xx[c, s][[2]], c], 
    bc1[xx[c, s][[3]], c], bc1[xx[c, s][[4]], c]}, {1, 1, 1, 
    1}, {bc4[xx[c, s][[1]], c], bc4[xx[c, s][[2]], c], 
    bc4[xx[c, s][[3]], c], bc4[xx[c, s][[4]], c]}, {bc2[xx[c, s][[1]]], 
    bc2[xx[c, s][[2]]], bc2[xx[c, s][[3]]], bc2[xx[c, s][[4]]]}};


(* Compute the determinant substituting $c$ and $s$ with explicit expressions.
   The expression for $c$ will depend on $s$ and lambda. The parameters = I omega is purely imaginary *)


det[lambda_, omega_, sp_] = Det[m[k + I omega e - (lambda  I omega)/(1 + I omega t),I omega, sp]];


(* Compute the solution )
solc = FindRoot[{Re[det[lambda, omega, oursp]] == 0, 
   Im[det[lambda, omega, oursp]] == 0}, {lambda, 
   7.2}, {omega, 290}]
(out: {lambda -> 7.22856, omega -> 285.74}*)

The problem of accuracy:

{lambdac,omegac} = {lambda,omega}/.solc;
det[lambdac, omegac, oursp]
(*out: -4.36557*10^-10 + 4.07454*10^-10 I*)

P.S. The scope of this computation is to find the non-trivial solutions A=(A1,A2,A3,A4)^T, of m(λc, ωc)A=0

P.P.S. I've tried to read the answer to Improve Accuracy of FindRoot, but I fear it is too complicated for my case.

Domen
  • 23,608
  • 1
  • 27
  • 45
fritess
  • 33
  • 2

1 Answers1

4

For this particular case, you can use exact numbers instead of inexact by applying Rationalize to all your parameters. Then, use Solve instead of NSolve. Finally, increase WorkingPrecision of FindRoot.

Clear[sp]
ll = 58.3 // Rationalize;
xin = 0.0034 // Rationalize;
b = 1730;
a = 0.185 // Rationalize;

...


root = Solve[eigen[c, s] == 0, x];


...


solc = FindRoot[{Re[det[lambda, omega, oursp]] == 0, 
   Im[det[lambda, omega, oursp]] == 0}, {lambda, 7.2}, {omega, 290}, 
  WorkingPrecision -> 30]


{lambdac, omegac} = {lambda, omega} /. solc;
det[lambdac, omegac, oursp]
(* {lambda -> 7.22856127683466811826154273590, 
    omega -> 285.740490749591546034984866881} )
( 0.10^-22 + 0.10^-22 I *)
Domen
  • 23,608
  • 1
  • 27
  • 45
  • Thanks! It solved the problem :) – fritess Feb 12 '24 at 15:47
  • 1
    @user123 Note that at machine precision, your original result of FindRoot was off from this solution by one or two units in the last place (a relative error from the true root that is less that 4 x 10^-16). It's pretty good and not really needing improvement, if you're going to use machine precision. The residual in this answer, less than 10^-22 starting with 30 digit precision, also shows that the original result, less than 10^-9 starting with 16 digit precision, has similar conditioning near the root and is about as good as can be expected. – Goofy Feb 12 '24 at 18:49