why the program is not running .

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

Answer 1

[Walter Roberson]

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

[neetu malik]

i dont understand sir,please correct my programe there is any mistake.

[Walter Roberson]

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]

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]

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]

Error using plot

Vectors must be the same length.

Error in g3 (line 78)

plot(k,D2,'b')

[Walter Roberson]

    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]

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]

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]

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.

[neetu malik]

this program is not stop. after run this program my leptop stop working.

[neetu malik]

please help in this program.

[Walter Roberson]

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]

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]

also if we disp(D) then the output is not come.

[neetu malik]

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]

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]

How can I solve this sir

[neetu malik]

How can I make graph between real(c) and k and imag(c) and k for different values of h

[Walter Roberson]

How can I solve this sir

Rent time on https://en.wikipedia.org/wiki/Frontier_(supercomputer) maybe.

[Walter Roberson]

I've been executing for 10 or more hours just to convert D to an anonymous function :(

[Walter Roberson]

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]

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]

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]

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

[neetu malik]

am I able to make graph between k and c by this program?

[neetu malik]

i want to graph between k and c for h=2.4, 2.6 and 2.8

[Walter Roberson]

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.