現在この質問をフォロー中です
- フォローしているコンテンツ フィードに更新が表示されます。
- コミュニケーション基本設定に応じて電子メールを受け取ることができます。
why the program is not running .
1 回表示 (過去 30 日間)
古いコメントを表示
tic
syms c h k
c11=106.8;
c33=54.57;
c13=9.68;
c15=0.28;
c35=-1.69;
c55=25.05;
r1=2727;
l=2.46;
u=5.66;
r2=7800;
i=sqrt(-1);
n=2.54;
a=(c33*c55);
b=(c35^2);
d=(c15*c33);
e=(c13*c35);
f=(c33+c55);
g=(c15*c35);
j=(c11*c33);
m=(c13^2);
o=(c13*c55);
p=(c55+c11);
q=((c11*c55)-c15^2);
r=(c15+c35);
t=(c11*c35-c13*c15);
a0=a-b;
a1=-2*i*(d-e);
a2=f*r1*(c^2)-2*g-j+m+2*o;
a3=-2*i*(r*r1*(c^2)-t);
a4=(r1^2)*(c^4)-p*r1*(c^2)+q;
s=[a0 a1 a2 a3 a4];
root=roots(s);
s1=root(1);
s2=root(2);
s3=root(3);
s4=root(4);
m1=-((c55)*(s1^2)-2*i*(c15)*(s1)+r1*(c^2)-(c11))/((c35)*(s1^2)-i*(c13+c55)*(s1)-(c15));
m2=-((c55)*(s2^2)-2*i*(c15)*(s2)+r1*(c^2)-(c11))/((c35)*(s2^2)-i*(c13+c55)*(s2)-(c15));
m3=-((c55)*(s3^2)-2*i*(c15)*(s3)+r1*(c^2)-(c11))/((c35)*(s3^2)-i*(c13+c55)*(s3)-(c15));
m4=-((c55)*(s4^2)-2*i*(c15)*(s4)+r1*(c^2)-(c11))/((c35)*(s4^2)-i*(c13+c55)*(s4)-(c15));
alpha=sqrt((l+2*u)/r2);
beta=sqrt(u/r2);
b0=(n^2)*((beta^2)/alpha^2);
b1=(n^2)*((1-((beta^2)/(alpha^2))^2))+n*(((c^2)/(alpha^2))-n)+n*((beta^2)/(alpha^2))*((c^2)/(beta^2)-n);
b2=((beta^2)/(alpha^2))*((c^2)/(alpha^2)-n)*((c^2)/(beta^2)-n);
p=[b0 0 b1 0 b2];
proot=roots(p);
p1=proot(1);
p2=proot(2);
p3=proot(3);
p4=proot(4);
n1=(((c^2)/(alpha^2))-n+n*((beta^2)/(alpha^2))*(p1^2))/(i*n*(1-((beta^2)/(alpha^2))*p1));
n2=(((c^2)/(alpha^2))-n+n*((beta^2)/(alpha^2))*(p2^2))/(i*n*(1-((beta^2)/(alpha^2))*p2));
K1=i*(c15)-m1*s1*(c35)+(i*m1-s1)*(c55);
K2=i*(c15)-m2*s2*(c35)+(i*m2-s2)*(c55);
K3=i*(c15)-m3*s3*(c35)+(i*m3-s3)*(c55);
K4=i*(c15)-m4*s4*(c35)+(i*m4-s4)*(c55);
K5=i*(c13)-m1*s1*(c33)+(i*m1-s1)*(c35);
K6=i*(c13)-m2*s2*(c33)+(i*m2-s2)*(c35);
K7=i*(c13)-m3*s3*(c33)+(i*m3-s3)*(c35);
K8=i*(c13)-m4*s4*(c33)+(i*m4-s4)*(c35);
e1=exp(k*s1*h);
e2=exp(k*s2*h);
e3=exp(k*s3*h);
e4=exp(k*s4*h);
A=[1 1 1 1 -1 -1; m1 m2 m3 m4 -n1 -n2; K1 K2 K3 K4 -n*u*(i*n1-p1) -n*u*(i*n2-p2); K5 K6 K7 K8 -n*(i*l-(l+2*u)*p1*n1) -n*(i*l-(l+2*u)*p2*n2); K1*e1 K2*e2 K3*e3 K4*e4 0 0; K5*e1 K6*e2 K7*e3 K8*e4 0 0];
AA =sym('AA',[6 6]);
detA=det(AA);
D=subs(detA, AA, A);
for h=2.4
k=4:0.2:5;
D1=solve(D==0,c);
D2=real(D1);
plot(k,D2,'b')
end
hold on
for h=2.6
k=4:0.2:5;
D1=solve(D==0,c);
D2=real(D1);
plot(k,D2,'r')
end
hold on
for h=2.8
k=4:0.2:5;
D1=solve(D==0,c);
D2=real(D1);
plot(k,D2,'r')
end
hold off
toc
採用された回答
Walter Roberson
2022 年 8 月 3 日
k=4:0.2:5;
D1=solve(D==0,c);
Consider the following code:
A = 1
B = 10*A
A = 2
What is the value of B after those three lines of code? Is B a formula that automatically updates as the original variables change? (If it is then B=B+1 leads to some interesting consequences)
When you syms k and assign an expression in k to a variable, and then assign a number to k, then does the variable get updated to reflect the new value of k? If so you can get into some contradictions easily.
You should read the documentation for subs()
27 件のコメント
Walter Roberson
2022 年 8 月 3 日
syms k
That is equivalent to
k = sym('k')
The symbolic engine creates a variable that lives inside the symbolic engine, and returns to matlab a reference to the symbol. MATLAB will receive something similar to '_ans[8859_42]' . MATLAB then stores that inside the MATLAB variable k
D = k
The reference _ans[8859_42] is what the matlab variable contains so the reference is copied to D
k = 4
The matlab level replaces the variable k with the numeric 4. The symbolic engine is not notified. The reference _ans etc is still valid and still refers to something inside the symbolic engine
D
D contains that reference _ans etc. The display routine sees that a symbolic expression is being used, and passes that _ans etc reference to the symbolic engine along with the instruction "return a text version of whatever this refers to". The symbolic engine sees that it refers to symbolic k and formats that k as text and returns 'k' to matlab to be displayed.
Notice here that matlab does not even know that D refers to k. There might have been a context such as k-k that the symbolic engine computed as not referring to k anymore. matlab does not know. matlab passes everything symbolic to the symbolic engine and stores references and asks the symbolic engine to do work, but matlab has no idea what the expressions mean. So when you replaced k at the MATLAB level with a number, matlab has no idea that might change symbolic expressions. MATLAB just has that '_ans[8859_42]' character string and passes it to the symbolic engine without knowing anything about it.
So when you assign k=4:.2:5 in your code, doing that has no effect on what D refers to. You are asking to solve() a version of D that still has the symbolic variable k inside it.
You need to read the documentation for the subs() call.
Walter Roberson
2022 年 8 月 3 日
編集済み: Walter Roberson
2022 年 8 月 3 日
You have another problem as well. If
syms R
k = 4:.2:5
D = R-k
solve(D==0,R)
then because k is a vector, R-k would be a vector, and D would be a vector involving R. The solve() call would then be asking to solve a vector of equations involving scalar R. But... solve() of a vector of equations asks solve() to find the value of the scalar variable R that satisfies all of the equations at the same time. It does not ask solve to find R for the first equation, and then find R for the second equation, and so on. solve() is for simultaneous equations.
This means that you have to loop over each of the numeric k values, asking to solve() for each one.
neetu malik
2022 年 8 月 3 日
編集済み: Walter Roberson
2022 年 8 月 4 日
tic
syms c h k
c11=106.8;
c33=54.57;
c13=9.68;
c15=0.28;
c35=-1.69;
c55=25.05;
r1=2727;
l=2.46;
u=5.66;
r2=7800;
i=sqrt(-1);
n=2.54;
a=sort((c33*c55));
b=sort((c35^2));
d=sort((c15*c33));
e=sort((c13*c35));
f=sort((c33+c55));
g=sort((c15*c35));
j=sort((c11*c33));
m=sort((c13^2));
o=sort((c13*c55));
p=sort((c55+c11));
q=sort(((c11*c55)-c15^2));
r=sort((c15+c35));
t=sort((c11*c35-c13*c15));
a0=a-b;
a1=-2*i*(d-e);
a2=f*r1*(c^2)-2*g-j+m+2*o;
a3=-2*i*(r*r1*(c^2)-t);
a4=(r1^2)*(c^4)-p*r1*(c^2)+q;
s=[a0 a1 a2 a3 a4];
root=roots(s);
s1=sort(root(1));
s2=sort(root(2));
s3=sort(root(3));
s4=sort(root(4));
m1=sort(-((c55)*(s1^2)-2*i*(c15)*(s1)+r1*(c^2)-(c11))/((c35)*(s1^2)-i*(c13+c55)*(s1)-(c15)));
m2=sort(-((c55)*(s2^2)-2*i*(c15)*(s2)+r1*(c^2)-(c11))/((c35)*(s2^2)-i*(c13+c55)*(s2)-(c15)));
m3=sort(-((c55)*(s3^2)-2*i*(c15)*(s3)+r1*(c^2)-(c11))/((c35)*(s3^2)-i*(c13+c55)*(s3)-(c15)));
m4=sort(-((c55)*(s4^2)-2*i*(c15)*(s4)+r1*(c^2)-(c11))/((c35)*(s4^2)-i*(c13+c55)*(s4)-(c15)));
alpha=sort(sqrt((l+2*u)/r2));
beta=sort(sqrt(u/r2));
b0=sort((n^2)*((beta^2)/alpha^2));
b1=sort((n^2)*((1-((beta^2)/(alpha^2))^2))+n*(((c^2)/(alpha^2))-n)+n*((beta^2)/(alpha^2))*((c^2)/(beta^2)-n));
b2=sort(((beta^2)/(alpha^2))*((c^2)/(alpha^2)-n)*((c^2)/(beta^2)-n));
p=[b0 0 b1 0 b2];
proot=roots(p);
p1=sort(proot(1));
p2=sort(proot(2));
p3=sort(proot(3));
p4=sort(proot(4));
n1=sort((((c^2)/(alpha^2))-n+n*((beta^2)/(alpha^2))*(p1^2))/(i*n*(1-((beta^2)/(alpha^2))*p1)));
n2=sort((((c^2)/(alpha^2))-n+n*((beta^2)/(alpha^2))*(p2^2))/(i*n*(1-((beta^2)/(alpha^2))*p2)));
K1=sort(i*(c15)-m1*s1*(c35)+(i*m1-s1)*(c55));
K2=sort(i*(c15)-m2*s2*(c35)+(i*m2-s2)*(c55));
K3=sort(i*(c15)-m3*s3*(c35)+(i*m3-s3)*(c55));
K4=sort(i*(c15)-m4*s4*(c35)+(i*m4-s4)*(c55));
K5=sort(i*(c13)-m1*s1*(c33)+(i*m1-s1)*(c35));
K6=sort(i*(c13)-m2*s2*(c33)+(i*m2-s2)*(c35));
K7=sort(i*(c13)-m3*s3*(c33)+(i*m3-s3)*(c35));
K8=sort(i*(c13)-m4*s4*(c33)+(i*m4-s4)*(c35));
e1=sort(exp(k*s1*h));
e2=sort(exp(k*s2*h));
e3=sort(exp(k*s3*h));
e4=sort(exp(k*s4*h));
A=[1 1 1 1 -1 -1; m1 m2 m3 m4 -n1 -n2; K1 K2 K3 K4 -n*u*(i*n1-p1) -n*u*(i*n2-p2); K5 K6 K7 K8 -n*(i*l-(l+2*u)*p1*n1) -n*(i*l-(l+2*u)*p2*n2); K1*e1 K2*e2 K3*e3 K4*e4 0 0; K5*e1 K6*e2 K7*e3 K8*e4 0 0];
AA =sym('AA',[6 6]);
detA=det(AA);
D=subs(detA, AA, A);
for h=2.4
syms R
k=4:0.2:5;
D=R-k;
D1=solve(D==0,R);
D2=real(D1);
plot(k,D2,'b')
end
hold on
for h=2.6
syms R
k=4:0.2:5;
D=R-k;
D1=solve(D==0,R);
D2=real(D1);
plot(k,D2,'g')
end
hold on
for h=2.8
syms R
k=4:0.2:5;
D=R-k;
D1=solve(D==0,R);
D2=real(D1);
plot(k,D2,'r')
end
hold off
toc
neetu malik
2022 年 8 月 3 日
Error using plot
Vectors must be the same length.
Error in g3 (line 78)
plot(k,D2,'b')
Walter Roberson
2022 年 8 月 4 日
syms R
R is a symbolic scalar
k=4:0.2:5;
k is a vector with 6 elements
D=R-k;
D will be a vector of 6 elements
D1=solve(D==0,R);
You are asking to solve the vector of 6 equations to find the single scalar value R that satisfies all of the equations simultaneously.
Walter Roberson
2022 年 8 月 4 日
syms x y
eqn1 = [3*x + 5*y == 7, 5*x + 7*y == 9]
eqn1 =
sol = solve(eqn1)
sol = struct with fields:
x: -1
y: 2
This asks to solve the vector of 2 equations. Why didn't it solve each one independently, like
for K = 1 : length(eqn1)
sols{K} = solve(eqn1(K), [x y]);
end
celldisp(sols)
sols{1} =
x: 7/3
y: 0
sols{2} =
x: 9/5
y: 0
Did you expect solve(eqn1) to solve each of the equations as distinct from each other, not interacting with the other? If not then why did you expect that
D1=solve(D==0,R);
would solve each element of D==0 independently of the other elements??
Walter Roberson
2022 年 8 月 4 日
編集済み: Walter Roberson
2022 年 8 月 4 日
tic
syms c h k
c11=106.8;
c33=54.57;
c13=9.68;
c15=0.28;
c35=-1.69;
c55=25.05;
r1=2727;
l=2.46;
u=5.66;
r2=7800;
i=sqrt(-1);
n=2.54;
a=sort((c33*c55));
b=sort((c35^2));
d=sort((c15*c33));
e=sort((c13*c35));
f=sort((c33+c55));
g=sort((c15*c35));
j=sort((c11*c33));
m=sort((c13^2));
o=sort((c13*c55));
p=sort((c55+c11));
q=sort(((c11*c55)-c15^2));
r=sort((c15+c35));
t=sort((c11*c35-c13*c15));
a0=a-b;
a1=-2*i*(d-e);
a2=f*r1*(c^2)-2*g-j+m+2*o;
a3=-2*i*(r*r1*(c^2)-t);
a4=(r1^2)*(c^4)-p*r1*(c^2)+q;
s=[a0 a1 a2 a3 a4];
root=roots(s);
s1=sort(root(1));
s2=sort(root(2));
s3=sort(root(3));
s4=sort(root(4));
m1=sort(-((c55)*(s1^2)-2*i*(c15)*(s1)+r1*(c^2)-(c11))/((c35)*(s1^2)-i*(c13+c55)*(s1)-(c15)));
m2=sort(-((c55)*(s2^2)-2*i*(c15)*(s2)+r1*(c^2)-(c11))/((c35)*(s2^2)-i*(c13+c55)*(s2)-(c15)));
m3=sort(-((c55)*(s3^2)-2*i*(c15)*(s3)+r1*(c^2)-(c11))/((c35)*(s3^2)-i*(c13+c55)*(s3)-(c15)));
m4=sort(-((c55)*(s4^2)-2*i*(c15)*(s4)+r1*(c^2)-(c11))/((c35)*(s4^2)-i*(c13+c55)*(s4)-(c15)));
alpha=sort(sqrt((l+2*u)/r2));
beta=sort(sqrt(u/r2));
b0=sort((n^2)*((beta^2)/alpha^2));
b1=sort((n^2)*((1-((beta^2)/(alpha^2))^2))+n*(((c^2)/(alpha^2))-n)+n*((beta^2)/(alpha^2))*((c^2)/(beta^2)-n));
b2=sort(((beta^2)/(alpha^2))*((c^2)/(alpha^2)-n)*((c^2)/(beta^2)-n));
p=[b0 0 b1 0 b2];
proot=roots(p);
p1=sort(proot(1));
p2=sort(proot(2));
p3=sort(proot(3));
p4=sort(proot(4));
n1=sort((((c^2)/(alpha^2))-n+n*((beta^2)/(alpha^2))*(p1^2))/(i*n*(1-((beta^2)/(alpha^2))*p1)));
n2=sort((((c^2)/(alpha^2))-n+n*((beta^2)/(alpha^2))*(p2^2))/(i*n*(1-((beta^2)/(alpha^2))*p2)));
K1=sort(i*(c15)-m1*s1*(c35)+(i*m1-s1)*(c55));
K2=sort(i*(c15)-m2*s2*(c35)+(i*m2-s2)*(c55));
K3=sort(i*(c15)-m3*s3*(c35)+(i*m3-s3)*(c55));
K4=sort(i*(c15)-m4*s4*(c35)+(i*m4-s4)*(c55));
K5=sort(i*(c13)-m1*s1*(c33)+(i*m1-s1)*(c35));
K6=sort(i*(c13)-m2*s2*(c33)+(i*m2-s2)*(c35));
K7=sort(i*(c13)-m3*s3*(c33)+(i*m3-s3)*(c35));
K8=sort(i*(c13)-m4*s4*(c33)+(i*m4-s4)*(c35));
e1=sort(exp(k*s1*h));
e2=sort(exp(k*s2*h));
e3=sort(exp(k*s3*h));
e4=sort(exp(k*s4*h));
A=[1 1 1 1 -1 -1; m1 m2 m3 m4 -n1 -n2; K1 K2 K3 K4 -n*u*(i*n1-p1) -n*u*(i*n2-p2); K5 K6 K7 K8 -n*(i*l-(l+2*u)*p1*n1) -n*(i*l-(l+2*u)*p2*n2); K1*e1 K2*e2 K3*e3 K4*e4 0 0; K5*e1 K6*e2 K7*e3 K8*e4 0 0];
AA =sym('AA',[6 6]);
detA=det(AA);
D=subs(detA, AA, A);
for h=2.4
syms R
k=4:0.2:5;
D=R-k;
D1 = arrayfun(@(X) solve(X==0,R), D);
D224=real(D1);
plot(k,D224,'b', 'displayname', '2.4')
end
hold on
for h=2.6
syms R
k=4:0.2:5;
D=R-k;
D1 = arrayfun(@(X) solve(X==0,R), D);
D226=real(D1);
plot(k,D226,'g', 'displayname', '2.6')
end
hold on
for h=2.8
syms R
k=4:0.2:5;
D=R-k;
D1 = arrayfun(@(X) solve(X==0,R), D);
D228=real(D1);
plot(k,D228,'r', 'displayname', '2.8')
end
hold off
legend show
toc
Elapsed time is 19.393399 seconds.
D224
D224 =
D226
D226 =
D228
D228 =
Notice that in the plot that you can only see a single line -- the last of the lines drawn. If you look at the D224, D226, D228 you can see that this is because the solutions to each of the solve() is exactly the same. This is because what you are asking to solve, R-k, is exactly the same for each h value. You do not ask to solve anything that involves h in the calculation, and you do not arrange that h gets substituted into any calculation to be solved.
Note that after you carefully calculate D according to det(), that you then overwrite D three times, with R-k each time.
neetu malik
2022 年 8 月 4 日
tic
syms c h k
c11=106.8;
c33=54.57;
c13=9.68;
c15=0.28;
c35=-1.69;
c55=25.05;
r1=2727;
l=2.46;
u=5.66;
r2=7800;
i=sqrt(-1);
n=2.54;
a=sort((c33*c55));
b=sort((c35^2));
d=sort((c15*c33));
e=sort((c13*c35));
f=sort((c33+c55));
g=sort((c15*c35));
j=sort((c11*c33));
m=sort((c13^2));
o=sort((c13*c55));
p=sort((c55+c11));
q=sort(((c11*c55)-c15^2));
r=sort((c15+c35));
t=sort((c11*c35-c13*c15));
a0=a-b;
a1=-2*i*(d-e);
a2=f*r1*(c^2)-2*g-j+m+2*o;
a3=-2*i*(r*r1*(c^2)-t);
a4=(r1^2)*(c^4)-p*r1*(c^2)+q;
s=[a0 a1 a2 a3 a4];
root=roots(s);
s1=sort(root(1));
s2=sort(root(2));
s3=sort(root(3));
s4=sort(root(4));
m1=sort(-((c55)*(s1^2)-2*i*(c15)*(s1)+r1*(c^2)-(c11))/((c35)*(s1^2)-i*(c13+c55)*(s1)-(c15)));
m2=sort(-((c55)*(s2^2)-2*i*(c15)*(s2)+r1*(c^2)-(c11))/((c35)*(s2^2)-i*(c13+c55)*(s2)-(c15)));
m3=sort(-((c55)*(s3^2)-2*i*(c15)*(s3)+r1*(c^2)-(c11))/((c35)*(s3^2)-i*(c13+c55)*(s3)-(c15)));
m4=sort(-((c55)*(s4^2)-2*i*(c15)*(s4)+r1*(c^2)-(c11))/((c35)*(s4^2)-i*(c13+c55)*(s4)-(c15)));
alpha=sort(sqrt((l+2*u)/r2));
beta=sort(sqrt(u/r2));
b0=sort((n^2)*((beta^2)/alpha^2));
b1=sort((n^2)*((1-((beta^2)/(alpha^2))^2))+n*(((c^2)/(alpha^2))-n)+n*((beta^2)/(alpha^2))*((c^2)/(beta^2)-n));
b2=sort(((beta^2)/(alpha^2))*((c^2)/(alpha^2)-n)*((c^2)/(beta^2)-n));
p=[b0 0 b1 0 b2];
proot=roots(p);
p1=sort(proot(1));
p2=sort(proot(2));
p3=sort(proot(3));
p4=sort(proot(4));
n1=sort((((c^2)/(alpha^2))-n+n*((beta^2)/(alpha^2))*(p1^2))/(i*n*(1-((beta^2)/(alpha^2))*p1)));
n2=sort((((c^2)/(alpha^2))-n+n*((beta^2)/(alpha^2))*(p2^2))/(i*n*(1-((beta^2)/(alpha^2))*p2)));
K1=sort(i*(c15)-m1*s1*(c35)+(i*m1-s1)*(c55));
K2=sort(i*(c15)-m2*s2*(c35)+(i*m2-s2)*(c55));
K3=sort(i*(c15)-m3*s3*(c35)+(i*m3-s3)*(c55));
K4=sort(i*(c15)-m4*s4*(c35)+(i*m4-s4)*(c55));
K5=sort(i*(c13)-m1*s1*(c33)+(i*m1-s1)*(c35));
K6=sort(i*(c13)-m2*s2*(c33)+(i*m2-s2)*(c35));
K7=sort(i*(c13)-m3*s3*(c33)+(i*m3-s3)*(c35));
K8=sort(i*(c13)-m4*s4*(c33)+(i*m4-s4)*(c35));
e1=sort(exp(k*s1*h));
e2=sort(exp(k*s2*h));
e3=sort(exp(k*s3*h));
e4=sort(exp(k*s4*h));
A=[1 1 1 1 -1 -1; m1 m2 m3 m4 -n1 -n2; K1 K2 K3 K4 -n*u*(i*n1-p1) -n*u*(i*n2-p2); K5 K6 K7 K8 -n*(i*l-(l+2*u)*p1*n1) -n*(i*l-(l+2*u)*p2*n2); K1*e1 K2*e2 K3*e3 K4*e4 0 0; K5*e1 K6*e2 K7*e3 K8*e4 0 0];
D=det(A);
for h=2.4
k=4:0.2:5;
D1 = solve(D==0,c);
D2=real(D1);
plot(k,D2,'b', 'displayname', '2.4')
end
hold on
for h=2.6
k=4:0.2:5;
D1 = solve(D==0,c);
D3=real(D1);
plot(k,D3,'g', 'displayname', '2.6')
end
hold on
for h=2.8
k=4:0.2:5;
D1 = solve(D==0,c);
D4=real(D1);
plot(k,D4,'r', 'displayname', '2.8')
end
hold off
legend show
i want to graph between c and k for differtent value of h. here D is the determinant of 6x6 matrix. when we take D=0 it become a equation in c,k,h. how can i solve this sir.
Walter Roberson
2022 年 8 月 4 日
D=det(A);
Taking the det() of a matrix of expressions is much slower than taking the det() of a representative matrix and subsituting in the coefficients later, the way I showed you with the AA matrix.
neetu malik
2022 年 8 月 4 日
編集済み: Walter Roberson
2022 年 8 月 4 日
tic
syms c h k
c11=106.8;
c33=54.57;
c13=9.68;
c15=0.28;
c35=-1.69;
c55=25.05;
r1=2727;
l=2.46;
u=5.66;
r2=7800;
i=sqrt(-1);
n=2.54;
a=sort((c33*c55));
b=sort((c35^2));
d=sort((c15*c33));
e=sort((c13*c35));
f=sort((c33+c55));
g=sort((c15*c35));
j=sort((c11*c33));
m=sort((c13^2));
o=sort((c13*c55));
p=sort((c55+c11));
q=sort(((c11*c55)-c15^2));
r=sort((c15+c35));
t=sort((c11*c35-c13*c15));
a0=a-b;
a1=-2*i*(d-e);
a2=f*r1*(c^2)-2*g-j+m+2*o;
a3=-2*i*(r*r1*(c^2)-t);
a4=(r1^2)*(c^4)-p*r1*(c^2)+q;
s=[a0 a1 a2 a3 a4];
root=roots(s);
s1=sort(root(1));
s2=sort(root(2));
s3=sort(root(3));
s4=sort(root(4));
m1=sort(-((c55)*(s1^2)-2*i*(c15)*(s1)+r1*(c^2)-(c11))/((c35)*(s1^2)-i*(c13+c55)*(s1)-(c15)));
m2=sort(-((c55)*(s2^2)-2*i*(c15)*(s2)+r1*(c^2)-(c11))/((c35)*(s2^2)-i*(c13+c55)*(s2)-(c15)));
m3=sort(-((c55)*(s3^2)-2*i*(c15)*(s3)+r1*(c^2)-(c11))/((c35)*(s3^2)-i*(c13+c55)*(s3)-(c15)));
m4=sort(-((c55)*(s4^2)-2*i*(c15)*(s4)+r1*(c^2)-(c11))/((c35)*(s4^2)-i*(c13+c55)*(s4)-(c15)));
alpha=sort(sqrt((l+2*u)/r2));
beta=sort(sqrt(u/r2));
b0=sort((n^2)*((beta^2)/alpha^2));
b1=sort((n^2)*((1-((beta^2)/(alpha^2))^2))+n*(((c^2)/(alpha^2))-n)+n*((beta^2)/(alpha^2))*((c^2)/(beta^2)-n));
b2=sort(((beta^2)/(alpha^2))*((c^2)/(alpha^2)-n)*((c^2)/(beta^2)-n));
p=[b0 0 b1 0 b2];
proot=roots(p);
p1=sort(proot(1));
p2=sort(proot(2));
p3=sort(proot(3));
p4=sort(proot(4));
n1=sort((((c^2)/(alpha^2))-n+n*((beta^2)/(alpha^2))*(p1^2))/(i*n*(1-((beta^2)/(alpha^2))*p1)));
n2=sort((((c^2)/(alpha^2))-n+n*((beta^2)/(alpha^2))*(p2^2))/(i*n*(1-((beta^2)/(alpha^2))*p2)));
K1=sort(i*(c15)-m1*s1*(c35)+(i*m1-s1)*(c55));
K2=sort(i*(c15)-m2*s2*(c35)+(i*m2-s2)*(c55));
K3=sort(i*(c15)-m3*s3*(c35)+(i*m3-s3)*(c55));
K4=sort(i*(c15)-m4*s4*(c35)+(i*m4-s4)*(c55));
K5=sort(i*(c13)-m1*s1*(c33)+(i*m1-s1)*(c35));
K6=sort(i*(c13)-m2*s2*(c33)+(i*m2-s2)*(c35));
K7=sort(i*(c13)-m3*s3*(c33)+(i*m3-s3)*(c35));
K8=sort(i*(c13)-m4*s4*(c33)+(i*m4-s4)*(c35));
e1=sort(exp(k*s1*h));
e2=sort(exp(k*s2*h));
e3=sort(exp(k*s3*h));
e4=sort(exp(k*s4*h));
A=[1 1 1 1 -1 -1; m1 m2 m3 m4 -n1 -n2; K1 K2 K3 K4 -n*u*(i*n1-p1) -n*u*(i*n2-p2); K5 K6 K7 K8 -n*(i*l-(l+2*u)*p1*n1) -n*(i*l-(l+2*u)*p2*n2); K1*e1 K2*e2 K3*e3 K4*e4 0 0; K5*e1 K6*e2 K7*e3 K8*e4 0 0];
AA =sym('AA',[6 6]);
detA=det(AA);
D=subs(detA, AA, A);
for h=2.4
k=4:0.2:5;
D1 = solve(D==0,c);
D2=real(D1);
plot(k,D2,'b', 'displayname', '2.4')
end
hold on
for h=2.6
k=4:0.2:5;
D1 = solve(D==0,c);
D3=real(D1);
plot(k,D3,'g', 'displayname', '2.6')
end
hold on
for h=2.8
k=4:0.2:5;
D1 = solve(D==0,c);
D4=real(D1);
plot(k,D4,'r', 'displayname', '2.8')
end
hold off
legend show
i also try this sir but yet this is not working
neetu malik
2022 年 8 月 4 日
Warning: Cannot find explicit solution.
> In solve (line 316)
In g1 (line 75)
Error using plot
Vectors must be the same length.
Error in g1 (line 77)
plot(k,D2,'b', 'displayname', '2.4')
this error is come, please help sir
Walter Roberson
2022 年 8 月 4 日
Finding a solution to this takes a number of hours. The expressions that are produced are quite large, and some of the sub-expressions are very finely balanced. The 288 sub-expressions making up the sum that is D each involve several exponentials multiplied together, with the exponentials each having a potential division by 0, and each of the sub-expressions involve complex numbers. To find a place where the solution is 0, you have to find a place where the real and imaginary components of all of the 288 parts all exactly balance out to 0. The solution probably involves letting the values of most of the sub-expressions more or less "float free" but carefully position c for the last sub-expression to be just on the edge of one of the discontinuities so that it can balance out the total contributions of the other 287 sub-expressions.
neetu malik
2022 年 8 月 4 日
How can I make graph between real(c) and k and imag(c) and k for different values of h
Walter Roberson
2022 年 8 月 5 日
I've been executing for 10 or more hours just to convert D to an anonymous function :(
Walter Roberson
2022 年 8 月 5 日
At about the 12 hours mark,
Error using symengine
Error: Expression is too large for MATLAB to evaluate.
so it cannot be done in one piece.
Walter Roberson
2022 年 8 月 5 日
It is feasible to use children() on D, and matlabFunction() on each of the parts, and then construct a function handle that invokes each of the parts and sums the result of the invocations.
However, invoking for even just one (complex) point takes 3 1/2 minutes...
Walter Roberson
2022 年 8 月 5 日
The good news (relatively) is that if you invoke the above reconstructed function on arrays of points, that although there is a lot of initial cost, that in bulk the timing is not terrible -- for example 700 seconds for 2500 points, about 1/4 second per point.
But you wouldn't want to do something like fzero() or fsolve(), as that invokes with individual points; ga 'usevectorized' has more potential for optimization, in theory.
Walter Roberson
2022 年 8 月 6 日
編集済み: Walter Roberson
2022 年 8 月 6 日
The symbolic approach is just not feasible. I've had it executing for over 10 hours just trying to find a good location for one particular h and k combination, and it is not close to a solution at all.
Using a numeric approach is much more satisfying.
This numeric code does not manage to solve the determinants in any of the tries. Instead, it gets the determinants down to about 3e-5 to 0.015 . At the moment we do not know if there are any actual solutions near the points found; it could be the case that there is an asymptope. They are areas worth investigating further.
Note that the fsolve() here often fail from any given starting point; when a failure is detected, I generate a random point nearby and try again. A decent point is usually found after a small number of tries.
T0 = tic;
syms c h k
c11=106.8;
c33=54.57;
c13=9.68;
c15=0.28;
c35=-1.69;
c55=25.05;
r1=2727;
l=2.46;
u=5.66;
r2=7800;
i=sqrt(-1);
n=2.54;
a=sort((c33*c55));
b=sort((c35^2));
d=sort((c15*c33));
e=sort((c13*c35));
f=sort((c33+c55));
g=sort((c15*c35));
j=sort((c11*c33));
m=sort((c13^2));
o=sort((c13*c55));
p=sort((c55+c11));
q=sort(((c11*c55)-c15^2));
r=sort((c15+c35));
t=sort((c11*c35-c13*c15));
a0=a-b;
a1=-2*i*(d-e);
a2=f*r1*(c^2)-2*g-j+m+2*o;
a3=-2*i*(r*r1*(c^2)-t);
a4=(r1^2)*(c^4)-p*r1*(c^2)+q;
s=[a0 a1 a2 a3 a4];
root=roots(s);
s1=sort(root(1));
s2=sort(root(2));
s3=sort(root(3));
s4=sort(root(4));
m1=sort(-((c55)*(s1^2)-2*i*(c15)*(s1)+r1*(c^2)-(c11))/((c35)*(s1^2)-i*(c13+c55)*(s1)-(c15)));
m2=sort(-((c55)*(s2^2)-2*i*(c15)*(s2)+r1*(c^2)-(c11))/((c35)*(s2^2)-i*(c13+c55)*(s2)-(c15)));
m3=sort(-((c55)*(s3^2)-2*i*(c15)*(s3)+r1*(c^2)-(c11))/((c35)*(s3^2)-i*(c13+c55)*(s3)-(c15)));
m4=sort(-((c55)*(s4^2)-2*i*(c15)*(s4)+r1*(c^2)-(c11))/((c35)*(s4^2)-i*(c13+c55)*(s4)-(c15)));
alpha=sort(sqrt((l+2*u)/r2));
beta=sort(sqrt(u/r2));
b0=sort((n^2)*((beta^2)/alpha^2));
b1=sort((n^2)*((1-((beta^2)/(alpha^2))^2))+n*(((c^2)/(alpha^2))-n)+n*((beta^2)/(alpha^2))*((c^2)/(beta^2)-n));
b2=sort(((beta^2)/(alpha^2))*((c^2)/(alpha^2)-n)*((c^2)/(beta^2)-n));
p=[b0 0 b1 0 b2];
proot=roots(p);
p1=sort(proot(1));
p2=sort(proot(2));
p3=sort(proot(3));
p4=sort(proot(4));
n1=sort((((c^2)/(alpha^2))-n+n*((beta^2)/(alpha^2))*(p1^2))/(i*n*(1-((beta^2)/(alpha^2))*p1)));
n2=sort((((c^2)/(alpha^2))-n+n*((beta^2)/(alpha^2))*(p2^2))/(i*n*(1-((beta^2)/(alpha^2))*p2)));
K1=sort(i*(c15)-m1*s1*(c35)+(i*m1-s1)*(c55));
K2=sort(i*(c15)-m2*s2*(c35)+(i*m2-s2)*(c55));
K3=sort(i*(c15)-m3*s3*(c35)+(i*m3-s3)*(c55));
K4=sort(i*(c15)-m4*s4*(c35)+(i*m4-s4)*(c55));
K5=sort(i*(c13)-m1*s1*(c33)+(i*m1-s1)*(c35));
K6=sort(i*(c13)-m2*s2*(c33)+(i*m2-s2)*(c35));
K7=sort(i*(c13)-m3*s3*(c33)+(i*m3-s3)*(c35));
K8=sort(i*(c13)-m4*s4*(c33)+(i*m4-s4)*(c35));
e1=sort(exp(k*s1*h));
e2=sort(exp(k*s2*h));
e3=sort(exp(k*s3*h));
e4=sort(exp(k*s4*h));
A=[1 1 1 1 -1 -1; m1 m2 m3 m4 -n1 -n2; K1 K2 K3 K4 -n*u*(i*n1-p1) -n*u*(i*n2-p2); K5 K6 K7 K8 -n*(i*l-(l+2*u)*p1*n1) -n*(i*l-(l+2*u)*p2*n2); K1*e1 K2*e2 K3*e3 K4*e4 0 0; K5*e1 K6*e2 K7*e3 K8*e4 0 0];
toc(T0)
T0 = tic;
Afun = matlabFunction(A, 'vars', [c, h, k]);
toc(T0)
Hvals = 2.4:0.2:2.8;
kvals = 4:0.2:5;
nH = length(Hvals);
nk = length(kvals);
results = zeros(nk, nH);
fvals = zeros(nk, nH);
opts = optimoptions('fsolve', 'Display', 'iter', 'MaxIter', 250, 'MaxFunctionEvaluations', 1000);
for hidx = 1 : nH
good = 0.296806441821817 - 0.00413848590710105i;
H = Hvals(hidx);
for kidx = 1 : nk
K = kvals(kidx);
Dfun = @(c) det(Afun(c, H, K));
errorflag = -inf;
while errorflag < 0
[results(kidx, hidx), fvals(kidx,hidx), errorflag] = fsolve(Dfun, good, opts);
if errorflag < 0
good = good + randn()/100+randn/1000*1i;
else
good = results(kidx, hidx);
fprintf('okay @(%d,%d)\n', kidx, hidx)
end
end
end
end
surf(Hvals, kvals, real(results))
Walter Roberson
2022 年 8 月 6 日
Yes. That surf() command at the bottom plots with H on the x axis, and k on the y axis, and c on the z axes.
Note that that particular surf() call only displays the real component of the c values. Since the values are complex, you will want to change how you choose to do the plotting.
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Get Started with Symbolic Math Toolbox についてさらに検索
タグ
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!エラーが発生しました
ページに変更が加えられたため、アクションを完了できません。ページを再度読み込み、更新された状態を確認してください。
Web サイトの選択
Web サイトを選択すると、翻訳されたコンテンツにアクセスし、地域のイベントやサービスを確認できます。現在の位置情報に基づき、次のサイトの選択を推奨します:
また、以下のリストから Web サイトを選択することもできます。
最適なサイトパフォーマンスの取得方法
中国のサイト (中国語または英語) を選択することで、最適なサイトパフォーマンスが得られます。その他の国の MathWorks のサイトは、お客様の地域からのアクセスが最適化されていません。
南北アメリカ
- América Latina (Español)
- Canada (English)
- United States (English)
ヨーロッパ
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
アジア太平洋地域
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)