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$
4
  • $\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
    Commented Jan 26, 2018 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
    Commented Jan 26, 2018 at 5:58
  • $\begingroup$ Sure FindRoot can give complex roots, e.g., FindRoot[x^2 - 3 I, {x, 1}]. $\endgroup$
    – Carl Woll
    Commented Jan 26, 2018 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
    Commented Jan 26, 2018 at 8:00

1 Answer 1

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$
5
  • $\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
    Commented Jan 26, 2018 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$ Commented Jan 26, 2018 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
    Commented Jan 27, 2018 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$ Commented Jan 27, 2018 at 15:01
  • $\begingroup$ @ Daniel Thanks $\endgroup$
    – T S Singh
    Commented Jan 29, 2018 at 16:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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