I am trying to do root-finding on a complex function. The function is a 3x3 matrix in the Laplacian domain. I would like to find the values of “s”, where the determinant of this matrix is zero. I am not sure what the best way to go about it is.
I will show the code, but first describe the problem a bit. More detail is found here particularly equation 12. When alpha (commented out in the code and substituted directly in Amat) is zero or close to zero, the determinant becomes a 5th order polynomial, which is easy to solve. In this case there is the zero frequency root, two more roots along with their complex conjugates for a total of 5 roots.
However when alpha becomes significant as when the fudge factor F is applied as per page 7 in the paper, then this is no longer the case. In fact according to the referenced paper there are an infinite number of roots because of the complex valued transcendental function exp(). What is a good algorithm to use to find the roots for this problem?
using Pkg
using LinearAlgebra
using Roots
"""
function DeterminantEqn12()
Polymer cover induced self-excited vibrations of nipped rolls
by Anssi T. Karttunen & Raimo von Hertzen
"""
# function DeterminantEqn12()
##
## Input simple spring-mass variables (no damping)
K = K_BottomFrame = 265.0e6 # N/m stiffness of frame or foundation at botttom
C = C_BottomRoll = 65400.0 # Ns/m damping of bottom roll
M = M_BottomRoll = 5380.0 # kg mass of bottom roll
k∞ = K_inf = 615.5e6 # N/m roll cover stiffess component
k₁ = K_1 = 99.28e6 # N/m roll cover stiffness next to dashpot
c₁ = C_1 = 84.38e3 # Ns/m SLS damping
τ = Tau_1 = C_1/K_1 # 0.85 ms for this example (4th line from bottom on pg 7)
m = M_TopRoll = 569.0 # kg mass of top roll
c = C_TopRoll = 8350 # Ns/m damping of top roll
k = K_TopFrame = 51.8e6 # N/m stiffness of frame at the top
F = 1.0/350 # Factor to account for slowly varying viscoelastic element pg 7
A = -F * K_inf/(K_inf+K_1)/Tau_1 # equation 4
##
# equation 11 for alpha
# kc = k_c = 1 # set for one revolution for relaxation of cover
# α = alpha = exp(kc*A*T)*exp(-kc*sᵢ*T) # s - Laplace transform variable; T - time for one revolution
# α = alpha = exp(A*T)*exp(-sᵢ*T) # with kc = 1
T = 0.1 # for 10 Hz rotation
sᵢ = 0+0im # Initial Value
Amat(s) = [m*s^2+(c+c₁)*s+k+k∞-k₁*exp(A*T)*exp(-s*T) -c₁*s+k₁*exp(A*T)*exp(-s*T) -k∞;
-c₁*s c₁*s+k₁ -k₁;
-k∞+k₁*exp(A*T)*exp(-s*T) -k₁*(1+exp(A*T)*exp(-s*T)) M*s^2+C*s+k₁+k∞+K]
f(s) = det(Amat(s))
# return f(sᵢ), Amat(sᵢ)
fzero(f, sᵢ)
find_zero(f, sᵢ)
# end