Inconsistency of the figure and value on the curve
3 views (last 30 days)
Show older comments
Hi all. I wrote a code to find the curve for trace and determinant. The problem is that when I zoom and choose one value on the curve and put it into the code, the value obtained for T is different from what is shown on the curve.
For example, here, for c=0.05826 on the curve is -4.06892e-07, while when code runs for this specific value, gives me T =-2.4350e-04.
Could you please let me know what the problem is that this happened and how I can solve it? Thanks in advance for any help
vplc=0.19;
delta=2.5;
tau_max=0.18;
Ktau=0.045;
kc=0.1;
kh=0.05;
Kbar=0.000015;
kp=0.15;
kb=0.4;
Kpm=0.15;
Ke=7;
vs=0.002;
ks=0.1;
kplc=0.055;
ki=2;
gamma=5.5;
kipr=0.18;
c=0.05826
%c =0:0.000001:0.3; %first I plotted the curves for this interval, zoom and choose a point on the curve and run the code for that specific vale
p=(vplc./ki).*(c.^2./((kplc).^2+c.^2));
h=kh.^4./(c.^4+kh.^4);
Po=(p.^2.*c.^4.*h)./(p.^2.*c.^4.*h.*(1+kb)+kb.*kp.^2.*(kc.^4+c.^4-(c.^4.*kh.^4./(c.^4+kh.^4))));
ct=(vs./(gamma.*kipr.*Po)).*(c.^2./(c.^2+ks.^2))+c.*(1+(1/gamma));
Podenominator=(p.^2.*c.^4.*h.*(1+kb)+kb.*kp.^2.*(kc.^4+c.^4-((c.^4.*kh.^4)./(c.^4+kh.^4))));
Dcnumerator=(8.*c.^7.*kh.^4*vplc.^2)./(ki.^2*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^2) - (4.*c.^9.*kh.^4.*vplc.^2)./(ki.^2.*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^3) - (4.*c.^11.*kh.^4*vplc.^2)./(ki.^2.*(c.^4 + kh.^4).^2.*(c.^2 + kplc.^2).^2);
Dcdenominator=kb.*kp.^2.*(4.*c.^3 - (4*c.^3*kh.^4)./(c.^4 + kh.^4) + (4.*c.^7*kh.^4)./(c.^4 + kh.^4).^2) + (8.*c.^7*kh.^4.*vplc.^2.*(kb + 1))./(ki.^2.*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^2) - (4*c.^9.*kh.^4.*vplc.^2.*(kb + 1))/(ki.^2.*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^3) - (4*c.^11.*kh^4.*vplc^2.*(kb + 1))./(ki^2.*(c.^4 + kh^4).^2.*(c.^2 + kplc.^2).^2);
DcPo=(Dcnumerator.*Podenominator-Dcdenominator.*(p.^2.*c.^4.*h))./(Podenominator.^2);
DhPo=((c.^4.*p.^2.*(Podenominator))-(c.^8.*p.^4.*(kh.^4./(c.^4+kh.^4)).*(1+kb)))./(Podenominator.^2);
Dcf1=kipr.*DcPo.*gamma.*ct-kipr.*DcPo.*(1+gamma).*c-kipr.*Po.*(gamma+1)-((2*vs.*c.*(ks^2))./(c.^2+ks^2).^2);
Dhf1=kipr.*DhPo.*(gamma.*ct-(1+gamma).*c);
Dcf2=(4.*c.^3./tau_max).*((kh.^8./(c.^4+kh.^4).^2))-(4.*c.^3./tau_max).*(kh.^4./(c.^4+kh.^4));
Dhf2=-c.^4./tau_max;
T = Dcf1 + Dhf2
D = Dcf1.*Dhf2-Dhf1.*Dcf2
eig1=(T+sqrt(T.^2-4.*D))./2;
eig2=(T-sqrt(T.^2-4.*D))./2;
figure
plot(c,T,'b')
hold on;
plot(c,D,'r')
hold on
plot(c,T.^2-4.*D,'g')
xlabel('c')
ylabel('T & D')
xlim([0 0.3])
ylim([-0.001 0.003])
0 Comments
Accepted Answer
Paul
on 18 Oct 2022
Edited: Paul
on 18 Oct 2022
Hi @M
The equation for Dcdenominator is not properly vectorized. In one spot it uses /, which should be ./
% was
% Dcdenominator=kb.*kp.^2.*(4.*c.^3 - (4*c.^3*kh.^4)./(c.^4 + kh.^4) + (4.*c.^7*kh.^4)./(c.^4 + kh.^4).^2) + (8.*c.^7*kh.^4.*vplc.^2.*(kb + 1))./(ki.^2.*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^2) - (4*c.^9.*kh.^4.*vplc.^2.*(kb + 1))/(ki.^2.*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^3) - (4*c.^11.*kh^4.*vplc^2.*(kb + 1))./(ki^2.*(c.^4 + kh^4).^2.*(c.^2 + kplc.^2).^2);
%
% should be
% xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx^use./
Dcdenominator=kb.*kp.^2.*(4.*c.^3 - (4*c.^3*kh.^4)./(c.^4 + kh.^4) + (4.*c.^7*kh.^4)./(c.^4 + kh.^4).^2) + (8.*c.^7*kh.^4.*vplc.^2.*(kb + 1))./(ki.^2.*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^2) - (4*c.^9.*kh.^4.*vplc.^2.*(kb + 1))./(ki.^2.*(c.^4 + kh.^4).*(c.^2 + kplc.^2).^3) - (4*c.^11.*kh^4.*vplc^2.*(kb + 1))./(ki^2.*(c.^4 + kh^4).^2.*(c.^2 + kplc.^2).^2);
More Answers (0)
See Also
Categories
Find more on Annotations 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!