Need help plotting this

Question

Amanda Lococo 2018-4-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/394754-need-help-plotting-this

评论： Star Strider 2018-4-15

I need help plotting this code so that my plot will look like the discontinuous sine waves (bolded black) in the graph posted below the code (for one unit cell). Right now my plot is showing up blank.

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];
B1 = [1 1;im*30 -im*30];
B2 = [1 1;im*5 -im*5];
C1a = [1 0;0 1];
C2a = [exp(im*k2*0.4) 0;0 exp(-im*k2*0.4)];
a1 = inv(B1)*V(:,1);
a2 = inv(C2a)*inv(B2)*T1*B1*a1;
for x1 = 0:0.1:0.4
    C1 = @(x1)([exp(im*k1*x1) 0;0 exp(-im*k1*x1)]);
    C1 = C1(x1);
    y1 = @(x1)(B1*C1*a1);
    y1 = y1(x1);
end
for x2 = 0.4:0.1:1
    C2 = @(x2)([exp(im*k2*x2) 0;0 exp(-im*k2*x2)]);
    C2 = C2(x2);
    y2 = @(x2)(B2*C2*a2);
    y2 = y2(x2);
end
plot(x1,real(y1(1,:)))
hold on
plot(x2,real(y2(1,:)))
xlabel('Position, x')
ylabel('Displacement, u')

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2018-4-13

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/394754-need-help-plotting-this#answer_314996

Two related problems are that the ‘y1’ and ‘y2’ functions do not use their arguments in their calculations. (I would also rename them ‘y1f’ and ‘y2f’ to avoid confusion with your ‘y1’ and ‘y2’ vectors.)

What do you want to do in those functions?

6 个评论
显示 4更早的评论隐藏 4更早的评论

Star Strider 2018-4-13

在 MATLAB Online 中打开

I have no idea what you’re doing. This is a bit more efficient, and produces plots:

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))/(1i*dTotal);
k2 = log(D(2,2))/(1i*dTotal);
k = [k1;k2];
B1 = [1 1;1i*30 -1i*30];
B2 = [1 1;1i*5 -1i*5];
C1a = [1 0;0 1];
C2a = [exp(1i*k2*0.4) 0;0 exp(-1i*k2*0.4)];
a1 = B1\V(:,1);
a2 = (B2*C2a)\T1*B1*a1;
x1 = 0:0.1:0.4;
C1 = @(x1)([exp(1i*k1*x1) 0;0 exp(-1i*k1*x1)]);
y1 = zeros(2, numel(x1));
for k = 1:numel(x1)
    y1(:,k) = (B1*C1(x1(k))*a1);
end
x2 = 0.4:0.1:1;
C2 = @(x2)([exp(1i*k2*x2) 0;0 exp(-1i*k2*x2)]);
y2 = zeros(2, numel(x2));
for k = 1:numel(x2)
    y2(:,k) = (B2*C2(x2(k))*a2);
end
plot(x1,real(y1(1,:)))
hold on
plot(x2,real(y2(1,:)))
xlabel('Position, x')
ylabel('Displacement, u')

I can’t completely vectorise this because I’m not certain what you’re doing. I moved the function definitions out of the loops, and called them in the ‘y1’ and ‘y2’ calculations.

You’ll have to determine if the plots are correct. I added preallocations for ‘y1’ and ‘y2’.