1
$\begingroup$

I am new to mathematica (less then a month), with the discussion available on this platform and with some tutorial available I have manage to write a code (not a smart one) to solve for the parameter values (values can be real or complex) for which the determinant is zero using Reduce function.

Objective: Search for value of 'sb' for which the det(DCoeffMat) is zero for different values of 'Pb'

(*%%% Parameters %%%%*)
ClearAll[sb, Pb, b, gama, L, w, h, Y, MI, EI, mu, rho, A, CN]
L = 10;
w = 0.02;
h = 0.02;
Y = 209*10^9;
MI = w*(h^3)/12;
EI = Y*MI;
mu = 0.001;
rho = 7800;
A = w*h;
CN = 4*Pi*mu/(Log[2*L/w] + 0.5) // N;
gama = rho*A*EI/(CN^2*L^4)
Do[B = b /. NSolve[b^4 + Pb*b^2 + sb + gama*sb^2 == 0, b] // N;
     DCoeffMat = {{1, 1, 1, 1}, {B[[1]]^3 + Pb*B[[1]], 
     B[[2]]^3 + Pb*B[[2]], B[[3]]^3 + Pb*B[[3]], 
     B[[4]]^3 + Pb*B[[4]]}, {B[[1]]*Exp[B[[1]]], B[[2]]*Exp[B[[2]]], 
     B[[3]]*Exp[B[[3]]], 
     B[[4]]*Exp[B[[4]]]}, {B[[1]]^2*Exp[B[[1]]] + Pb*Exp[B[[1]]], 
     B[[2]]^2*Exp[B[[2]]] + Pb*Exp[B[[2]]], 
     B[[3]]^2*Exp[B[[3]]] + Pb*Exp[B[[3]]], 
     B[[4]]^2*Exp[B[[4]]] + Pb*Exp[B[[4]]]}};
 DA = Det[DCoeffMat];
 Print[Reduce[
     DA == 0 && -1000 < Re[sb] < 1000 && -1000 < Im[sb] < 1000, sb] //
     Quiet // N]
 , {Pb, 0.1, 4., 0.1}]

Here is the equation used for finding B Polynomial for solving **"b"**

Here is the matrix whose determinant is equated to zero enter image description here

On running the program its showing an error.

Am I missing something???

As per the literature from where the problem is pick up it says that for pb<2.4674 the roots of sb has small real part (with non zero imaginary part) and after that its real value increases.

Thanks in advance.

$\endgroup$
  • $\begingroup$ You can't use NSolve with sb as a parameter. You're missing a ; after the // N and before DCoeffMat. And it's unlikely that Reduce will be able to solve your equation. You will most likely need to use FindRoot instead. $\endgroup$ – Carl Woll Jan 26 '18 at 5:49
  • $\begingroup$ Thanks for the correction. I found that NSolve can give the four roots of the polynomial equation. Can Findroot function give complex roots??? $\endgroup$ – T S Singh Jan 26 '18 at 5:58
  • $\begingroup$ Sure FindRoot can give complex roots, e.g., FindRoot[x^2 - 3 I, {x, 1}]. $\endgroup$ – Carl Woll Jan 26 '18 at 6:25
  • $\begingroup$ @Carl, I have replaced the solver with the code FindRoot[DA, {sb, I}] , its working. But I want vary the initial point (in the above case I have given 0+0*I). How can I do that within the Do loop?? With the motive to find different solutions. $\endgroup$ – T S Singh Jan 26 '18 at 8:00
2
$\begingroup$

Not a good idea to start variable names with capitals, hence the renaming below.

ll = 10;
w = 0.02;
h = 0.02;
yy = 209*10^9;
mi = w*(h^3)/12;
ei = yy*mi;
mu = 0.001;
rho = 7800;
aa = w*h;
cn = 4*Pi*mu/(Log[2*ll/w] + 0.5) // N;
gama = rho*aa*ei/(cn^2*ll^4);

Now define the functions of interest so they only operate on explicit numeric input.

bB[pb_?NumericQ, sb_?NumericQ] := 
  b /. NSolve[b^4 + pb*b^2 + sb + gama*sb^2 == 0, b];
dCoeffMat[pb_?NumericQ, sb_?NumericQ] := 
 With[{b = bB[pb, sb]}, {{1, 1, 1, 1}, {b[[1]]^3 + pb*b[[1]], 
    b[[2]]^3 + pb*b[[2]], b[[3]]^3 + pb*b[[3]], 
    b[[4]]^3 + pb*b[[4]]}, {b[[1]]*Exp[b[[1]]], b[[2]]*Exp[b[[2]]], 
    b[[3]]*Exp[b[[3]]], 
    b[[4]]*Exp[b[[4]]]}, {b[[1]]^2*Exp[b[[1]]] + pb*Exp[b[[1]]], 
    b[[2]]^2*Exp[b[[2]]] + pb*Exp[b[[2]]], 
    b[[3]]^2*Exp[b[[3]]] + pb*Exp[b[[3]]], 
    b[[4]]^2*Exp[b[[4]]] + pb*Exp[b[[4]]]}}]
dA[pb_?NumericQ, sb_?NumericQ] := Det[dCoeffMat[pb, sb]]

From here one can use FindRoot to get a value for sb that makes the determinat vanish given a value for pb. I use the real part of the determinant to avoid imaginary fuzz.

sbvals = 
 sb /. Table[FindRoot[Re[DA[pb, sb]] == 0, {sb, 1}], {pb, .1, 4., .1}]

(* Out[4217]= {0.0000893249467205, 0.000180282234353, 0.000271244538972, \
0.00036220809792, 0.000453172158726, 0.000544136470395, \
0.000635100925702, 0.000726065470127, 0.000817030074505, \
0.000907994720715, 0.000998959397339, 0.00108992409677, \
0.00118088881389, 0.00127185354457, 0.00136281828633, \
0.00145378303708, 0.0015447477952, 0.00163571255949, \
0.00172667732895, 0.00181764210281, 0.00190860688046, \
0.00199957166138, 0.00209053644524, 0.00218150123126, \
0.00227246601965, 0.00236343081001, 0.002454395602, 0.00254536039557, \
0.00281298528292, 0.00312767097318, 0.00341574060368, \
0.00368348997729, 0.00393511426993, 0.00417357535177, \
0.00440105694776, 0.00461922434613, 0.00482938236888, \
0.00503257643349, 0.00522965991584, 0.00542134060187} *)
$\endgroup$
  • $\begingroup$ @ Daniel Thanks. The results shows real values of sb. I am interested in both the Re[sb] and Im[sb]. How and where do I have to modify in your code ? $\endgroup$ – T S Singh Jan 26 '18 at 18:22
  • $\begingroup$ Try removing the Re[...] wrapper in that FindRoot and maybe give initial guess {sb, 1+I} (I don't have time to test this at the moment). $\endgroup$ – Daniel Lichtblau Jan 26 '18 at 18:31
  • $\begingroup$ @ Daniel I have implemented your suggestion, now the imaginary parts of sb are visible. One last doubt (hope not troubling you much :-)) suppose in implementing FindRoot if I want to vary the initial guess (as complex number) with the purpose to search roots in wide range of complex plane. How can it be done?? $\endgroup$ – T S Singh Jan 27 '18 at 3:22
  • $\begingroup$ You could use a Table that accounts for all the starting values. Something like Table[FindRoot[Re[DA[pb, sb]] == 0, {sb, sbval}], {pb, .1, 4., .1}], {sbval, sbinits}] where sbinits is the list of initial values you want to try. $\endgroup$ – Daniel Lichtblau Jan 27 '18 at 15:01
  • $\begingroup$ @ Daniel Thanks $\endgroup$ – T S Singh Jan 29 '18 at 16:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.