# Select positive value from a vector

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]]


• Use something like Pick[soln, Reduce[ForAll[{k, m}, k > 0 && m > 0, Positive[#]]] & /@ soln] (adjust restrictions as needed).
– ciao
May 17, 2014 at 23:51
• @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. May 18, 2014 at 1:56
• 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.
– ciao
May 18, 2014 at 2:06
• @Kuba Sure, I'm going to post it. May 24, 2014 at 17:43

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}]