Inconsistency of the figure and value on the curve

2 次查看（过去 30 天）
M
M 2022-10-17
评论： Paul 2022-10-18
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])

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

采纳的回答

Paul
Paul 2022-10-18
编辑：Paul 2022-10-18
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);
  2 个评论
M
M 2022-10-18
Excellent @Paul, it works perfectly... thanks
I have a question, I would appreciate if you can answer it.
With another method that I trust, I found a point where T~0 (when c~0.091), but even though all my calculations (DcPo, DhPo, T and D...)and parameters are correct here (I checked them in Maple), it contradicts with what I got (for example here for c~0,058, T becomes almost zero). you don't have any suggestion how to find the problem?

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Exploration and Visualization 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by