Need Help Plotting Mode Shapes
18 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Prajit T R
2018-4-9
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 个评论
Prajit T R
2018-4-10
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.
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!