Need Help Plotting Mode Shapes
13 views (last 30 days)
Show older comments
I need this code to plot mode shapes, but my plots are coming up blank. Thanks in advance!
clear all;
format long;
im = sqrt(-1);
CellLength = 1;
ibeta = 1;
%Define materal properties
CellLength = 1;
layers = 2;
d = [0.4;0.6];
dTotal = d(1,1)+d(2,1);
xc = [0;0.4];
Ef = 12;
pf = 3;
cf = sqrt(Ef/pf);
Em = 1;
pm = 1;
cm = sqrt(Em/pm);
w = 5;
T1 = [cos(.2*w) (1/(6*w))*sin(.2*w); -6*w*sin(.2*w) cos(.2*w)];
T2 = [cos(.6*w) (1/w)*sin(.6*w); -w*sin(.6*w) cos(.6*w)];
T = T2*T1;
Z1 = 6*w;
Z2 = w;
Z = [Z1;Z2];
%Solve eigenvalue problem for k
[V,D] = eig(T); %D = eigenvalues, %V = eigenvectors
k1 = log(D(1,1))/(im*dTotal);
k2 = log(D(2,2))/(im*dTotal);
k = [k1;k2];
for j = 1:layers
B = @(j)([1 1; im*Z(j,1) -im*Z(j,1)]);
B = B(j);
C_a = @(j)([exp(im*k(j,1)*xc(j,1)) 0; 0 exp(-im*k(j,1)*xc(j,1))]);
C_a = C_a(j);
if j == 1
a = (inv(B)*V(:,1));
alpha = a;
beta = B;
else
a = inv(C_a)*inv(B)*T*beta*alpha;
end
for x = 0:0.1:5
C = @(x)([exp(im*k(j,1)*x) 0; 0 exp(-im*k(j,1)*x)]);
C = C(x);
y = @(x)(B*C*a);
y = y(x);
end
end
plot(x,real(y(:,1)))
0 Comments
Accepted Answer
Prajit T R
on 9 Apr 2018
Hi Amanda
I am assuming that you wish to plot the variation of y against x for the values from 0 to 5 in steps of 0.5 as per the loop above. This can be done using the following code:
clear all;
format long;
im = sqrt(-1);
CellLength = 1;
ibeta = 1;
%Define materal properties
CellLength = 1;
layers = 2;
d = [0.4;0.6];
dTotal = d(1,1)+d(2,1);
xc = [0;0.4];
Ef = 12;
pf = 3;
cf = sqrt(Ef/pf);
Em = 1;
pm = 1;
cm = sqrt(Em/pm);
w = 5;
T1 = [cos(.2*w) (1/(6*w))*sin(.2*w); -6*w*sin(.2*w) cos(.2*w)];
T2 = [cos(.6*w) (1/w)*sin(.6*w); -w*sin(.6*w) cos(.6*w)];
T = T2*T1
Z1 = 6*w;
Z2 = w;
Z = [Z1;Z2];
%Solve eigenvalue problem for k
[V,D] = eig(T); %D = eigenvalues, %V = eigenvectors
k1 = log(D(1,1))/(im*dTotal);
k2 = log(D(2,2))/(im*dTotal);
k = [k1;k2];
for j = 1:layers
B = @(j)([1 1; im*Z(j,1) -im*Z(j,1)]);
B = B(j);
C_a = @(j)([exp(im*k(j,1)*xc(j,1)) 0; 0 exp(-im*k(j,1)*xc(j,1))]);
C_a = C_a(j);
if j == 1
a = (inv(B)*V(:,1));
alpha = a;
beta = B;
else
a = inv(C_a)*inv(B)*T*beta*alpha;
end
L=[];
M=[];
for x = 0:0.1:5
C = @(x)([exp(im*k(j,1)*x) 0; 0 exp(-im*k(j,1)*x)]);
C = C(x);
y = @(x)(B*C*a);
y = y(x);
L(end+1)=y(1);
M(end+1)=y(2);
end
end
plot(0:0.1:5,real(L))
hold on
plot(0:0.1:5,real(M))
hold off
In this code, L is the equivalent for y(:,1) that you were trying to achieve. M on the other hand is y(:,2). I hope this helps.
Cheers
2 Comments
Prajit T R
on 10 Apr 2018
y=y(x) inside the loop actually contains two elements. You can observe that by removing the semicolon next to the statement. So, what L(end+1) does is to append the first element of y to a list L, which is initially empty. Similarly, the second element of y is appended to M.
At the end of the iterations, L and M would contain all the complex values generated.
More Answers (0)
See Also
Categories
Find more on Line Plots in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!