Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
2  
Use something like Pick[soln, Reduce[ForAll[{k, m}, k > 0 && m > 0, Positive[#]]] & /@ soln] (adjust restrictions as needed). –  rasher May 17 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. –  Lawerance May 18 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. –  rasher May 18 at 2:06
    
@Kuba Sure, I'm going to post it. –  Lawerance May 24 at 17:43
add comment

1 Answer 1

up vote 1 down vote accepted

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}]
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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