-
Notifications
You must be signed in to change notification settings - Fork 1
/
string_eigenvalue_solver.m
67 lines (52 loc) · 1.8 KB
/
string_eigenvalue_solver.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
function [ naturalFrequencies, naturalModes ] = ...
string_eigenvalue_solver( masses, F, h )
% string_eigenvalue_solver Solves simple eigenvalue problem of a string
% Solves simple eigenvalue problem of a string with user-given masses in
% a column or row vector, constant tension force F on the string and
% separation of masses h
% Get dimension of arrays by size of masses vector
% Getting the maximum of the dimensions will allow the mass vector
% to be a column vector or row vector
n = max(size(masses,1), size(masses,2));
% --- Generate A ---
% 2s on main diagonal
A = 2*eye(n);
% -1s on second diagonals
diagOnes = -1 * ones(n-1,1);
A = A + diag(diagOnes,1) + diag(diagOnes,-1);
% --- Generate D ---
% Diagonal elements
diags = F ./ (h * masses);
D = diag(diags);
% Get eigenvalues and eigenvectors (natrual frequencies) of D * A
[naturalModes, eigenVals] = eig(D * A);
% Gather naturalModes in case a gpuArray was used
naturalModes = gather(naturalModes);
% Natural frequencies are square roots of eigenvalues times -i
naturalFrequencies = -1i * sqrt(eigenVals);
% --- Plotting ---
% Range of numbers for plotting
maxVal = h * n;
x = h:h:maxVal;
% Create new figure
plotFig = figure;
plotAx = axes;
hold(plotAx);
% Markers for plotting
markers = {'--+',':o','-.*','--x',':s','-.d', ...
'--^',':v','-.>','--<',':p','-.h'};
% Plot eigenmodes with frequencies
for i=1:size(naturalModes, 2)
plot(plotAx, x, naturalModes(:,i), markers{mod(i,numel(markers)) + 1}, ...
'DisplayName', sprintf('Mode %d, Frequency = %f', ...
i, eigenVals(i,i)), 'Linewidth', 2);
end
% Set limits of plot
xlim(plotAx, [min(x) - 1, max(x) + 1]);
% Show legend
plotLegend = legend('-DynamicLegend');
set(plotLegend, 'FontSize', 30);
set(plotLegend, 'Location', 'southoutside');
% Toggle off hold on axes
hold(plotAx);
end