how can i make graph between real(c) and k for diffferent values of h
1 view (last 30 days)
Show older comments
tic
format shortG
syms c 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;
n=2.54;
i=sqrt(-1);
c1=c^2;
a01=c33*c55;
a02=c35^2;
a0=a01-a02;
a11=-2*c15*c33;
a12=2*c13*c35;
a1=(a11+a12)*i;
a111=c33+c55;
a211=(a111*r1);
a212=2*c15*c35;
a213=c11*c33;
a214=c13^2;
a215=2*c13*c55;
a21=vpa(2.1712e+05)*c1;
a22=-a212-a213+a214+a215;
a2=vpa(a21-5248.5);
a31=-2*(c15+c35)*r1;
a32=2*c11*c35;
a33=-2*c13*c15;
a3=vpa(((7690.1)*c1-360.98-5.4208)*i);
a41=r1^2;
a42=-(c55+c11)*r1;
a43=c11*c55;
a44=-c15^2;
a4=vpa(7436529*c1^2-(3.5955e+05)*c1+ 2675.3-0.0784);
S=[a0 a1 a2 a3 a4];
R=roots(S);
s11=R(1);
s1=vpa(s11);
s12=R(2);
s2=vpa(s12);
s13=R(3);
s3=vpa(s13);
s14=R(4);
s4=vpa(s14);
m11=c55*(s1^2)-2*i*c15*s1+r1*(c^2)-c11;
m12=c35*(s1^2)-i*(c13+c55)*s1-c15;
m1=-(m11/m12);
m21=c55*(s2^2)-2*i*c15*s2+r1*(c^2)-c11;
m22=c35*(s2^2)-i*(c13+c55)*s2-c15;
m2=-(m21/m22);
m31=c55*(s3^2)-2*i*c15*s3+r1*(c^2)-c11;
m32=c35*(s3^2)-i*(c13+c55)*s3-c15;
m3=-(m31/m32);
m41=c55*(s4^2)-2*i*c15*s4+r1*(c^2)-c11;
m42=c35*(s4^2)-i*(c13+c55)*s4-c15;
m4=-(m41/m42);
alpha1=(l+2*u)/r2;
alpha=vpa(alpha1,4);
beta1=u/r2;
beta=vpa(beta1,4);
d11=(n^2)*(beta/alpha);
d1=vpa(d11,4);
d21=(n^2)*((1-(beta/alpha))^2)+n*(((c^2)/alpha)-n)+n*(beta/alpha)*(((c^2)/beta)-n);
d2=vpa(d21,4);
d31=(beta/alpha)*(((c^2)/alpha)-n)*(((c^2)/beta)-n);
d3=vpa(d31,4);
P=[d1 0 d2 0 d3];
R1=roots(P);
p11=R1(1);
p1=vpa(p11,4);
p12=R1(2);
p2=vpa(p12,4);
n11=((c^2)/alpha)-n+n*(beta/alpha)*(p1^2);
n12=i*n*(1-(beta/alpha))*p1;
n1=n11/n12;
n21=((c^2)/alpha)-n+n*(beta/alpha)*(p2^2);
n22=i*n*(1-(beta/alpha))*p2;
n2=n21/n22;
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=-n*u*(i*n1-p1);
k6=-n*u*(i*n2-p2);
k7=i*c13-m1*s1*c33+(i*m1-s1)*c35;
k8=i*c13-m2*s2*c33+(i*m2-s2)*c35;
k9=i*c13-m3*s3*c33+(i*m3-s3)*c35;
k10=i*c13-m4*s4*c33+(i*m4-s4)*c35;
k11=-n*(i*l-(l+2*u)*p1*n1);
k12=-n*(i*l-(l+2*u)*p2*n2);
for h=2.4
e1=exp(k*s1*h);
e2=exp(k*s2*h);
e3=exp(k*s3*h);
e4=exp(k*s4*h);
%A=[m1+n2 m2+n2 m3+n2 m4+n2 n2-n1; k1+k6 k2+k6 k3+k6 k4+k6 k5-k6; k7+k12 k8+k12 k9+k12 k10+k12 k11-k12; k1*e1 k2*e2 k3*e3 k4*e4 0; k7*e1 k8*e2 k9*e3 k10*e4 0];
A1=[k1+k6 k2+k6 k3+k6 k4+k6; k7+k12 k8+k12 k9+k12 k10+k12; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
A2=[m1+n2 m2+n2 m3+n2 m4+n2; k7+k12 k8+k12 k9+k12 k10+k12; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
A3=[m1+n2 m2+n2 m3+n2 m4+n2; k1+k6 k2+k6 k3+k6 k4+k6; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
A1A1 =sym('A1A1',[4 4]);
detA1=det(A1A1);
D1=subs(detA1, A1A1, A1);
A2A2 =sym('A2A2',[4 4]);
detA2=det(A2A2);
D2=subs(detA2, A2A2, A2);
A3A3 =sym('A3A3',[4 4]);
detA3=det(A3A3);
D3=subs(detA3, A3A3, A3);
D11=(n2-n1)*D1;
D21=(k5-k6)*D2;
D31=(k11-k12)*D3;
D=D11-D21+D31;
for k=4:0.2:5
Dfun=matlabFunction(D);
result=solve(Dfun==0,c);
plot(k,real(result))
end
end
hold on
for h=2.6
e1=exp(k*s1*h);
e2=exp(k*s2*h);
e3=exp(k*s3*h);
e4=exp(k*s4*h);
%A=[m1+n2 m2+n2 m3+n2 m4+n2 n2-n1; k1+k6 k2+k6 k3+k6 k4+k6 k5-k6; k7+k12 k8+k12 k9+k12 k10+k12 k11-k12; k1*e1 k2*e2 k3*e3 k4*e4 0; k7*e1 k8*e2 k9*e3 k10*e4 0];
A1=[k1+k6 k2+k6 k3+k6 k4+k6; k7+k12 k8+k12 k9+k12 k10+k12; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
A2=[m1+n2 m2+n2 m3+n2 m4+n2; k7+k12 k8+k12 k9+k12 k10+k12; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
A3=[m1+n2 m2+n2 m3+n2 m4+n2; k1+k6 k2+k6 k3+k6 k4+k6; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
A1A1 =sym('A1A1',[4 4]);
detA1=det(A1A1);
D1=subs(detA1, A1A1, A1);
A2A2 =sym('A2A2',[4 4]);
detA2=det(A2A2);
D2=subs(detA2, A2A2, A2);
A3A3 =sym('A3A3',[4 4]);
detA3=det(A3A3);
D3=subs(detA3, A3A3, A3);
D11=(n2-n1)*D1;
D21=(k5-k6)*D2;
D31=(k11-k12)*D3;
D=D11-D21+D31;
for k=4:0.2:5
Dfun=matlabFunction(D);
result=solve(Dfun==0,c);
plot(k,real(result))
end
end
hold on
for h=2.8
e1=exp(k*s1*h);
e2=exp(k*s2*h);
e3=exp(k*s3*h);
e4=exp(k*s4*h);
%A=[m1+n2 m2+n2 m3+n2 m4+n2 n2-n1; k1+k6 k2+k6 k3+k6 k4+k6 k5-k6; k7+k12 k8+k12 k9+k12 k10+k12 k11-k12; k1*e1 k2*e2 k3*e3 k4*e4 0; k7*e1 k8*e2 k9*e3 k10*e4 0];
A1=[k1+k6 k2+k6 k3+k6 k4+k6; k7+k12 k8+k12 k9+k12 k10+k12; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
A2=[m1+n2 m2+n2 m3+n2 m4+n2; k7+k12 k8+k12 k9+k12 k10+k12; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
A3=[m1+n2 m2+n2 m3+n2 m4+n2; k1+k6 k2+k6 k3+k6 k4+k6; k1*e1 k2*e2 k3*e3 k4*e4; k7*e1 k8*e2 k9*e3 k10*e4];
A1A1 =sym('A1A1',[4 4]);
detA1=det(A1A1);
D1=subs(detA1, A1A1, A1);
A2A2 =sym('A2A2',[4 4]);
detA2=det(A2A2);
D2=subs(detA2, A2A2, A2);
A3A3 =sym('A3A3',[4 4]);
detA3=det(A3A3);
D3=subs(detA3, A3A3, A3);
D11=(n2-n1)*D1;
D21=(k5-k6)*D2;
D31=(k11-k12)*D3;
D=D11-D21+D31;
for k=4:0.2:5
Dfun=matlabFunction(D);
result=solve(Dfun==0,c);
plot(k,real(result))
end
end
hold off
toc
i try to draw the graph between real© and k for different value of h. but my programe is not run. please help
0 Comments
Answers (3)
Walter Roberson
on 16 Aug 2022
Replace
Dfun=matlabFunction(D);
with
Dfun = subs(D);
in multiple places.
2 Comments
Walter Roberson
on 16 Aug 2022
Oh yes, that code. It is not possible to get that code to work.
I posted numeric code that does find solutions, several days ago, at https://www.mathworks.com/matlabcentral/answers/1773150-why-the-program-is-not-running#comment_2304010 .
There is NO chance that solve() will return a symbolic solution for you.
vpasolve() might potentially return a solution for you, if it were given a guess very close to the actual solution. But even given a very good starting point, it would take hours to find one solution.
Your current approach is not practical. The numeric approach, on the other hand, is practical.
neetu malik
on 16 Aug 2022
1 Comment
Walter Roberson
on 16 Aug 2022
What is your numeric implementation of the current Question, rewritten along the lines of the numeric code that I posted in https://www.mathworks.com/matlabcentral/answers/1773150-why-the-program-is-not-running#comment_2304010 ?
neetu malik
on 17 Aug 2022
2 Comments
Walter Roberson
on 17 Aug 2022
Okay, go ahead and do that. But my advice from having worked with your equations is that you have no chance of doing this using solve() and no realistic chance of doing this using vpasolve() either. However if you write your code in numeric form the way I showed in https://www.mathworks.com/matlabcentral/answers/1773150-why-the-program-is-not-running#comment_2304010 then you have a chance.
The code I posted there already plots c against h and k in a 3D plot.
See Also
Categories
Find more on Directed Graphs in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!