-1
$\begingroup$

I have written a small program to find the natural frequencies of some system using an energy approach. v1 and t1 are the potentials and kinetic energy of the beam. v2 and t2 are the potentials and kinetic energy of the bar. v2 ad t2 depends upon γ, and γ is varying. So I have created a module to extract the t2 and v2 for every value of γ. and v3=v2 and t3=t2. and I have constructed lagrangian Lg=((t1+t2+t3)-(v1+v2+v3))+other terms. for every value of γ, we have an Lg now, this Lg is minimized WRT some coefficients and construct a matrix out ou it. For this, I have written a small module called sol. Now my program is not running and gives an error "The solution set contains a full-dimensional component". How to fix this?

 Y = 2*^11;
ρ = 7850;
aa = 0.1*0.1;
Iyy = 0.1^4/12;
L1 = 4;
z[1] = 1.5;
z[2] = 3;


v1 = 833333. (0.25367 c[1]^2 + 0.388654 c[1] c[2] + 1.21287 c[2]^2 + 
     0.297836 c[1] c[3] - 2.67104 c[2] c[3] + 2.17733 c[3]^2 + 
     0.0215263 c[1] c[4] - 0.622222 c[2] c[4] + 2.08561 c[3] c[4] + 
     15.2699 c[4]^2);
t1 = 39.25 ω^2 (0.666667 c[1]^2 + 1.02142 c[1] c[2] + 
     0.557505 c[2]^2 + 0.782739 c[1] c[3] + 0.129908 c[2] c[3] + 
     0.562914 c[3]^2 + 0.056483 c[1] c[4] - 0.286262 c[2] c[4] + 
     0.539437 c[3] c[4] + 0.664979 c[4]^2);
sd1 = Subdivide[0.001, 0.5, 300];
sd2 = Subdivide[0.55, 2, 300];
sd3 = Subdivide[2.2, 5, 50];
sd=Flatten[{sd1,sd2,sd3}];


γ = L1*sd;
fixedfree1 = 
  Table[Cos[(i*π*x2)/γ[[i]]], {i, 1, Length[γ]}];
fixedfree2 = 
  Table[Cos[(i*π*x2)/γ[[i]]], {i, 1, Length[γ]}];

U1 = Table[(Subscript[d, 1][1]*fixedfree1[[i]] + 
     Subscript[d, 1][2]*fixedfree2[[i]]), {i, 1, Length[γ]}];
U1x = Table[D[(U1[[i]]), {x2, 1}], {i, 1, Length[U1]}];

in3 = Table[Expand[(U1x[[i]])^2], {i, 1, Length[U1x]}];
in4 = Table[Expand[(U1[[i]])^2], {i, 1, Length[U1x]}];

var2 = {Subscript[d, 1][1], Subscript[d, 1][2]};


rulesPE = 
  Table[CoefficientRules[in3[[i]], var2], {i, 1, Length[γ]}];
g1[j_] := 
 Module[{potenergy}, rules3 = rulesPE[[j]]; 
  rules3[[All, 2]] = 
   Quiet[NIntegrate[rules3[[All, 2]], {x2, 0, γ[[j]]}]] // 
    Chop; potenergy = 0.5*Y*aa*(FromCoefficientRules[rules3, var2]); 
  Return[potenergy]]
v2 = Table[g1[i], {i, 1, Length[γ]}];

rulesKE = Table[CoefficientRules[in4[[i]], var2], {i, 1, Length[in4]}];
g2[j_] := 
 Module[{kinenergy}, rules4 = rulesKE[[j]]; 
  rules4[[All, 2]] = 
   Quiet[NIntegrate[rules4[[All, 2]], {x2, 0, γ[[j]]}]] // 
    Chop; kinenergy = 
   0.5*ρ*aa*ω^2*(FromCoefficientRules[rules4, var2]); 
  Return[kinenergy]]
t2 = Table[g2[i], {i, 1, Length[γ]}];

U2 = U1 /. {Subscript[d, 1] -> Subscript[d, 2], x2 -> x3};
v3 = v2 /. Subscript[d, 1] -> Subscript[d, 2];
t3 = t2 /. Subscript[d, 1] -> Subscript[d, 2];
var3 = var2 /. Subscript[d, 1] -> Subscript[d, 2];
var1 = Table[c[i], {i, 1, 4}]

T = Table[Simplify[t1 + t2[[i]] + t3[[i]]], {i, 1, Length[t2]}];
V = Table[Simplify[v1 + v2[[i]] + v3[[i]]], {i, 1, Length[v2]}];

W = Total[Table[c[i] Sin[(i*π*x1)/L1], {i, 1, 4}]];

Lg = Table[((T[[i]] - 
       V[[i]]) + λ1*((W /. x1 -> z[1]) - (U1[[i]] /. 
          x2 -> γ[[i]])) + λ2*((W /. 
          x1 -> z[2]) - (U2[[i]] /. x3 -> γ[[i]]))), {i, 1, 
    Length[γ]}];
Dimensions[Lg]

variables = Flatten[{var1, var2, var3, λ1, λ2}]
sol[j_] := 
  Module[{root}, 
   equations = 
    Table[D[Lg[[j]], {variables[[i]], 1}], {i, 1, Length[variables]}];
    Rarz = Normal@CoefficientArrays[equations, variables][[2]]; 
   P = Chop[Det[Rarz]]; s1 = NSolve[P == 0 && 0 < ω < 2000]; 
   s2 = ω /. s1; s3 = Surd[(ρ*aa*s2^2*L1^4)/(Y*Iyy), 4]; 
   s4 = s2/(2 π); root = s4; Return[root]];
nonbeta = Table[sol[i], {i, 1, Length[γ]}];
$\endgroup$
  • 1
    $\begingroup$ A modicum of debugging would indicate that those determinants evaluate to 0.0, and NSolve then gives that result: In[1707]:= NSolve[0 == 0 && 0 < x < 2000] During evaluation of In[1707]:= NSolve::fulldim: The solution set contains a full-dimensional component; use Reduce for complete solution information. Out[1707]= {{}}. I'm not voting to close this because working around the issue is nontrivial. But I will reiterate that some debugging would have lead to a healthier question a la "How to redress this zero determinant problem?" $\endgroup$ – Daniel Lichtblau Sep 20 at 14:46
1
$\begingroup$

The code is too long for me to try to go over all of it, encapsulate variables properly, and the like. The main issue is that Det for a matrix containing both symbolic and approximate numeric entries is problematic, and leads to results of 0.0 in cases I have tried. A way to work around this is to interpolate the determinants, using roots of unity as values for the variable. The code below seems to work well for this purpose.

myDet[mat_?MatrixQ] /; Length[Variables[mat]] == 1 := 
 myDet[mat, First[Variables[mat]]]

myDet[mat_?MatrixQ, x_] /; Variables[mat] === {x} := Catch[Module[
   {dims = Dimensions[mat], xvals, dvals, ipoly},
   If[! Apply[Equal, dims], Throw[$Failed]];
   xvals = Exp[2*I*Pi*Range[0, dims[[1]]]/(dims[[1]] + 1)];
   dvals = Table[Det[mat /. x -> v], {v, xvals}];
   ipoly = InterpolatingPolynomial[Transpose[{xvals, dvals}], x];
   Expand[ipoly/(ipoly /. x -> 0)]
   ]]

--- edit ---

Oops. I posted too quickly and right now do not have time to repair. Here is the issue: I did not get a proper upper bound on the degree of the determinant. A correction would be to find such a bound and use it in place of dims[[1]] where the interpolation values are produced.

--- end edit ---

$\endgroup$

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.