0 个投票

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')

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

 采纳的回答

Star Strider
Star Strider 2018-4-13

1 个投票

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更早的评论

Amanda Lococo
Amanda Lococo 2018-4-13
I wasn't sure if I needed the functions for those or not. I figured I didn't since the only x-iteration is in the function C2. I am trying to write a code based on the results of a paper and the function described is that which is stated: y(x) = BjCj(x)aj to be specific.
What I need the function to do overall is plot x against y, where y is changing as a function of x (because C depends on x).
Star Strider
Star Strider 2018-4-13
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’.
Amanda Lococo
Amanda Lococo 2018-4-15
The plots do work, but the plots I need are posted in the original thread. It should resemble a sinusoidal wave with discontinuities at the layer boundaries (in this case at 0.4). I appreciate all of your help!
Star Strider
Star Strider 2018-4-15
My pleasure.
If you extend ‘x2’ to at least 3, you get the sine curve:
x2 = 0.4:0.1:3;
and with ‘x1’ defined as it is, the discontinuity appears at 0.4.
Now that your code produces plots, you need to experiment with it to get the results you want. I have no idea what you are doing, so I can only help with the code, not the concept (that may be outside my areas of expertise).
Amanda Lococo
Amanda Lococo 2018-4-15
I will try that then! Thanks again for all of your help!
Star Strider
Star Strider 2018-4-15
As always, my pleasure!

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by