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 is there a way to use loop, like Do, For, to rewrite the last line in the code. I tried to use Do to make this work, but mathematica keeps telling me that "svars will be suppressed during this calculation". I want to do the loop way, not the Table way i used below, because Tabel gives me a really ugly format of eigenvector.

enter image description here enter image description here enter image description here

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]

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

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

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

eigenval = Pick[soln, Resolve[ForAll[{k, m}, k > 0 && m > 0, Positive[#]]] & /@ soln];

Do[Print["\n\nThe positive eigenfrequencies we computed are " \[Omega][i] , " = ", eigenval[[i]]], {i, 4}]

b = Table[a[i], {i, Nmax}];

Table[Solve[matrix.b == 0 /. \[Omega] -> eigenval[[i]], b], {i, 1,Nmax}]
share|improve this question
    
You mean something like Table[Solve[matrix.b == 0, b], {ω, eigenval}] /. Rule -> Equal // MatrixForm ? –  xzczd May 19 at 7:07
    
@xzczd Yes, but to use the Loop to find the eigenvectors. –  Lawerance May 19 at 7:11
    
Why "loop"? I think the output is not "ugly" anymore, and Table is essentially a loop, too. Or you prefer this format?: b /. Table[Solve[matrix.b == 0, b], {ω, eigenval}] // MatrixForm –  xzczd May 19 at 7:14
    
Do you know why it won't work whenI change Table to Do? –  Lawerance May 19 at 7:25
    
I'm not sure what you mean: Do[Solve[matrix.b == 0, b] /. Rule -> Equal // Transpose // MatrixForm // Print // Quiet, {ω, eigenval}]; Do[b /. Solve[matrix.b == 0, b] // Transpose // Quiet // MatrixForm // Print, {ω, eigenval}] –  xzczd May 19 at 7:32

1 Answer 1

up vote 2 down vote accepted

According to the added information for the original problem, I think your problem can be simply solved by:

Clear[k, m]
coe = k/m SparseArray[{Band[{1, 1}] -> -2, Band[{1, 2}] -> 1, Band[{2, 1}] -> 1}, {4, 4}];
efre = Sqrt[-coe // Eigenvalues] // Simplify
evec = coe // Eigenvectors // Simplify

If you insist on making the result "look nicer":

"The eigenfrequencies are:"
Column[Subscript[ω, #] & /@ Range@4 == efre // Thread, Spacings -> 2]

enter image description here

"The eigenvectors are:"
Column[Subscript[v, #] & /@ Range@4 == MatrixForm /@ evec // Thread, Spacings -> 2]

enter image description here

share|improve this answer
    
Thanks a lot! Great help from you! –  Lawerance May 19 at 9:30
    
@Lawerance Glad you like it :) –  xzczd May 19 at 9:34
    
I took a deep look at your code this afternoon, and came up with several questions. First, it seems Center does not have any effect on the alignment of eigenvectors, so why do you put it above? Second, can you explain this to me: MatrixForm /@ evec // Thread. So, you Map the MatrixForm to every evec. But when I try Map[evec] in the code, I thought it would be the same but I did not work, why? –  Lawerance May 19 at 20:29
    
@Lawerance Center does have no effect here, I just forgot to remove it 囧, edited. For the second question, Map[evec] doesn't make any sense, the FullForm of MatrixForm /@ evec is Map[MatrixForm, evec]. also, notice what's Threaded here is == i.e. the precedence is Thread[((Subscript[v, #1] & ) /@ Range[4]) == (MatrixForm /@ evec)], you can use Shift+Ctrl+. to check the precedence. –  xzczd May 20 at 3:00
    
Thanks! Just got a chance to look at your comments after finishing my midterm yesterday. –  Lawerance May 23 at 19:56

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.