I am using an modified version of Wagon's FindAllCrossings2D[] function available here Updating Wagon's FindAllCrossings2D[] function. I am defining f[alpha,w] and g[alpha,w] as the real and imaginary part of function da[alpha,w] whose complex roots are to be evaluated.
(*Code for Roots Finding*)
(*Options[FindAllCrossings2D] =
Sort[Join[
Options[FindRoot], {MaxRecursion -> Automatic,
PerformanceGoal :> $PerformanceGoal, PlotPoints -> Automatic}]];
FindAllCrossings2D[funcs_, {x_, xmin_, xmax_}, {y_, ymin_, ymax_},
opts___] :=
Module[{contourData, seeds, tt,
fy = Compile[{x, y}, Evaluate[funcs[[2]]]]},
contourData =
Map[First,
Cases[Normal[
ContourPlot[funcs[[1]], {x, xmin, xmax}, {y, ymin, ymax},
Contours -> {0}, ContourShading -> False,
PlotRange -> {Full, Full, Automatic},
Evaluate[
Sequence @@
FilterRules[Join[{opts}, Options[FindAllCrossings2D]],
DeleteCases[Options[ContourPlot], Method -> _]]]]], _Line,
Infinity]];
seeds =
Flatten[Map[#[[
1 + Flatten[
Position[Rest[tt = Sign[Apply[fy, #, 2]]] Most[tt], -1]]]] &,
contourData], 1];
If[seeds == {}, seeds,
Select[Union[
Map[{x, y} /.
FindRoot[{funcs[[1]] == 0,
funcs[[2]] == 0}, {x, #[[1]]}, {y, #[[2]]},
Evaluate[
Sequence @@
FilterRules[Join[{opts}, Options[FindAllCrossings2D]],
Options[FindRoot]]]] &,
seeds]], (xmin < #[[1]] < xmax && ymin < #[[2]] < ymax) &]]]*)
(*Program For Sliding Control*)
Clear[sigma, sigmab, xib, cn, ct, a, l, ei, tau, pb, fb, fc, ks, aa,bb, k, lamda, xi, b, lp, w, alpha, da, an]
cn = 0.0034;
ct = 0.0017;
a = (0.185*10^(-6))/2;
l = 58.3*10^(-6);
ei = 1700*10^(-24) // N;
tau = 0.004;
pb = 0.03;
fb = 3.8*10^(-12) // N;
fc = 2*10^(-12) // N;
fp = 1.8*10^(-6) // N;
ks = 94.8*10^(-3) // N;
gama = 0.273*10^(-3) // N;
rho = 150*10^(6) // N;
k3 = 800*10^(12) // N;
aa = 2*rho*pb*(1 - pb)*fb*fp/fc;
bb = 2*rho*pb*fp;
k = -1620;
lamda = -7.60;
sigma[alpha_, w_] := alpha + I*w;
xi[alpha_,
w_] := (-aa*(alpha + I*w)/(1 + (alpha + I*w)*tau)) + (bb*(alpha +
I*w));
sigmab[alpha_, w_] := (alpha + I*w)*cn*l^4/ei ;
xib[alpha_,
w_] := ((-aa*(alpha + I*w)/(1 + (alpha + I*w)*tau)) + (bb*(alpha +
I*w)))*a^2*l^2/ei;
lp[alpha_,
w_] := ((-aa*(alpha + I*w)/(1 + (alpha + I*w)*tau)) + (bb*(alpha +
I*w)))*a^2*l^2/ei;
nomodes = 1;
(*Calculation of Elements of Matrix*)
b = {1.875, 4.694, 7.8548, 10.9955, 14.1372, 17.2788, 20.4204,
23.5619, 26.7035, 29.8451, 32.9867, 36.1283, 39.2699, 42.4115,
45.5531, 48.6947};
an = (-Sin[b] + Sinh[b])/(Cos[b] + Cosh[b]);
q = an*(Cos[b*x] - Cosh[b*x]) + (Sin[b*x] + Sinh[b*x]);
qx = D[q, x];
qxx = D[qx, x];
qxxx = D[qxx, x];
phi = an*(-Cos[b*x] - Cosh[b*x]) + (-Sin[b*x] + Sinh[b*x]);
phix = D[phi, x];
phixx = D[phix, x];
phixxx = D[phixx, x];
cm = Table[
NIntegrate[phi[[i]]*q[[j]], {x, 0, 1}], {i, nomodes}, {j, nomodes}];
coeffMat[alpha_, w_] := Table[
k1m = NIntegrate[phixx[[i]]*qxx[[j]], {x, 0, 1}];
k2m = NIntegrate[phix[[i]]*qx[[j]], {x, 0, 1}];
k3ma = phi[[i]]*qxxx[[j]] + phi[[i]]*lp[alpha, w]*qx[[j]] -
phix[[i]]*qxx[[j]] /. x -> 1;
k3mb = phi[[i]]*qxxx[[j]] + phi[[i]]*lp[alpha, w]*qx[[j]] -
phix[[i]]*qxx[[j]] /. x -> 0;
k3m = k3ma - k3mb;
km = k1m + lp[alpha, w]*k2m + k3m, {i, nomodes}, {j, nomodes}]
da[alpha_, w_] = Det[coeffMat[alpha, w] + sigmab[alpha, w]*cm]
f[alpha_, w_] = Re[da[alpha, w]]
g[alpha_, w_] = Im[da[alpha, w]]
pts = FindAllCrossings2D[{f[alpha, w], g[alpha, w]}, {alpha, -100,
100}, {w, -400, 400},
Method -> {"Newton", "StepControl" -> "LineSearch"},
PlotPoints -> 85, WorkingPrecision -> 20] // Chop
ContourPlot[{f[alpha, w], g[alpha, w]}, {alpha, -100, 100}, {w, -400,
400}, Contours -> {0}, ContourShading -> False,
Epilog -> {AbsolutePointSize[6], Red, Point /@ pts}]
There are certain problems which I am facing
- While calculating matrix form of
coeffMat[alpha,w]usingTablefunction its showing error$Failed When using the
FindAllCrossings2D[]it is unable to find the roots(I think the first error is leading to the failure in finding the roots because of improper definition of the functions
coeffMat[alpha,w])
Please help me out!
Thanks in advance
f[alpha, w]cannot evaluateIm[w]orRe[alpha](and others) because it did not specify thatwandalphaareReals. Compare yourf[alpha, w]withFullSimplify[ComplexExpand@Re[f[alpha, w]], alpha > 0 && w > 0]for example. – anderstood Feb 06 '18 at 15:29