0
$\begingroup$

I am wondering how to pick out the positive Omega value from soln, and make these positive values into a row vector. I know how to do it if elements in soln are numerical. But here all the elements are in the symbolic form.

Remove["Global`*"]

Nmax = 4;

T = 1/2 Sum[m[i] x[i]'[t]^2, {i, 1, Nmax}];
U = 1/2 Sum[k[i] (x[i][t] - x[i - 1][t])^2, {i, 1, Nmax + 1}];
L = T - U;

EL[q_] := D[L, q] - D[D[L, D[q, t]], t]

eigen = Table[EL[x[i][t]], {i, 1, Nmax}];

x[i_][t_] = a[i] E^(I \[Omega] t);

sim = eigen /. t -> 0;
x[0][t] = 0;
x[Nmax + 1][t] = 0;

For[i = 1, i < Nmax + 1, i++, m[i] = m]
For[i = 1, i < Nmax + 2, i++, k[i] = k]

matrix = D[sim, {Array[a, Nmax]}];

Print["The reduced eigenmatrix looks like ", matrix // MatrixForm]

soln = \[Omega] /. Solve[Det[matrix] == 0, \[Omega]]

enter image description here

$\endgroup$
4
  • 2
    $\begingroup$ Use something like Pick[soln, Reduce[ForAll[{k, m}, k > 0 && m > 0, Positive[#]]] & /@ soln] (adjust restrictions as needed). $\endgroup$
    – ciao
    May 17, 2014 at 23:51
  • $\begingroup$ @rasher I feel a little confused here about Reduce. I look it up and it says that Reduce is used to solve and simplify the expression, basically like Solve. I feel like here Resolve rather than Reduce makes much more sense to me. $\endgroup$
    – Lawerance
    May 18, 2014 at 1:56
  • $\begingroup$ Equivalent for this use - Resolve attempts to remove quantifiers, not always solving. Reduce automagically applies Resolve as step 1. Use either here - they're equivalent for fully quantified systems. $\endgroup$
    – ciao
    May 18, 2014 at 2:06
  • $\begingroup$ @Kuba Sure, I'm going to post it. $\endgroup$
    – Lawerance
    May 24, 2014 at 17:43

1 Answer 1

1
$\begingroup$

Thanks for the help of @rasher and @xzczd, I'm now pretty much finishing this problem. I didn't use the Pick and Reduce command in this problem ultimately, because I found we can manually select the positive eigenvalues by taking the positive square root. But I still appreciate for the help given by these two people. Here is my code.

Remove["Global`*"]

Nmax = 4;

T = 1/2 Sum[m[i] x[i]'[t]^2, {i, 1, Nmax}];
U = 1/2 Sum[k[i] (x[i][t] - x[i - 1][t])^2, {i, 1, Nmax + 1}];
L = T - U;

EL[q_] := D[L, q] - D[D[L, D[q, t]], t]

eigen = Table[EL[x[i][t]], {i, 1, Nmax}];

x[i_][t_] = a[i] E^(I \[Omega] t);

sim = eigen /. t -> 0;
x[0][t] = 0;
x[Nmax + 1][t] = 0;

For[i = 1, i <= Nmax + 1, i++, m[i] = m; k[i] = k; m[5] = 0]

(* In order to factor out the vector a, we use the Jacobian \
transformation *)

matrix = D[sim, {Array[a, Nmax]}];

"Problem 3,4."
Print["\nThe reduced eigenmatrix looks like ", matrix // MatrixForm]

(* By looking at the reduced matrix, we got the intuitive sense that \
this matrix is symmetric; and more importantly, it's with dignoal and \
sub-diagonal elements *)

(* So, to simplify the code, we can actually construct the original \
matrix and then apply `Eigensystem` to get the eigenvalues, \
eigenvectors *)

org = k/m SparseArray[{Band[{1, 1}] -> 2, Band[{1, 2}] -> -1, 
     Band[{2, 1}] -> -1}, {4, 4}];
efre = Sqrt[
   org // Eigenvalues]; (* This step is crucial for sorting out the \
positive eigenvalues, therefore reducing the expected "8" to "4" \
eigenvalues *)
evec = org // Eigenvectors;

"\n\nThe positive eigenfrequencies are:\n\n"
Column[Subscript[\[Omega], #] & /@ Range @4 == efre // Thread, 
 Spacings -> 2]

"\n\nThe corresponding eigenvectors are:\n\n"
Column[Subscript[a, #] & /@ Range@4 == MatrixForm /@ evec // Thread, 
 Spacings -> 2]

m = 1.;
k = 1.;

sep = 2.;  
ysep = 2.;

range = {{-sep, sep*(Nmax + 1)}, {-ysep, (Nmax + 1)*ysep}};
tlo = 0;
thi = 100;

xsoln[t_] = 
  Table[sep*(mass - 1) + 
    evec[[mode, mass]] Cos[efre[[mode]] t], {mode, 1, Nmax}, {mass, 1,
     Nmax}];

plist[t_] = 
  Table[{PointSize[0.02], 
    Point[{xsoln[t][[mode, mass]], (mode - 1)*ysep}]}, {mode, 1, 
    Nmax}, {mass, 1, Nmax}];
plist[t_] = Flatten[plist[t], 2];   
modeplot[t_] := Graphics[plist[t]];

pix[t_] := Show[modeplot[t], Axes -> True, PlotRange -> range];

Animate[pix[t], {t, tlo, thi}]
$\endgroup$

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.