Need Help Plotting Mode Shapes

21 次查看(过去 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)))

采纳的回答

Prajit T R
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 个评论
Amanda Lococo
Amanda Lococo 2018-4-9
编辑:Amanda Lococo 2018-4-9
Thank you for this! Can you explain these two lines to me?
L(end+1)=y(1);
M(end+1)=y(2);
Thanks again!
Prajit T R
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!

Translated by