方程解取正实解问题。

大侠们帮我看看我这个解方程取正实解的问题呀，我看了好多网上的帖子取正值的，可我的取不了呀，我觉得是我的解是符号型号的问题，可是不会调呀，请大侠们帮调一下，代码如下：

syms p1 p2 t0 t2 u K L fai2 bai xm zz z

%%%%%%%%%%%%%%%%%%%%%%

K=1.4;u=(K-1)/(K+1);

%%%%%%%%% P2点 %%%%%%%

px2=3;py2=3;pz2=0;Q=1;

%%%%%%%%%%%%%%%%%%%%%%%% 爆心S坐标

sx2=9;sy2=3;sz2=1.75;p0=0.1013;ca=340;

%%%%%%%%%%%%%%%%%%%%%%%%%%%

R2=sqrt((sx2-px2)^2+(sy2-py2)^2+(sz2-pz2)^2) %%%%%%%%%%%%%%% 入射距离

Z2=R2/Q^(1/3)

%%%%%%%%%%%%%%%%%%%%% 入射角

fai2=acos(sz2/R2)*180/pi

%%%%%%%%%%%%%%%%%%%

p_i2=0.084/Z2+0.27/Z2^2+0.7/Z2^3

pai1_2=(p_i2+p0)/p0;

%%%%%%%%%%%%%%%%简化

Mi=sqrt(p_i2/(p0*(1+u))+1)

bai_JX2=1.75/(Mi-1)+39

%%%%%%%%%%%%%%%%简化

%%%%%%%%%%%%%%计算值

M=(pai1_2-1)*t0/(pai1_2+u+(1+pai1_2*u)*t0^2);

arfa=-M^2*(1-u)^2+2*M*(1-u)*t0+1;

beta=-M^2*(1-u)^2*t0-M*(1-2*u-t0^2)+t0;

tuga=M*(1+t0^2);

y=arfa^2-4*beta*tuga;

%%y=eval(x^3-Bc*x^2/Ac-Cc*x/Ac-Dc/Ac) %%赋值语句

x=solve(y)

vpa(x,6) %%%%

n=length(x);

%xc=[];

for j=1:numel(x)

if isreal(x(j))

x(j);

end

end

x

vpa(xc,6) %%%%

原来用下面这个就可以取正了，可是发现当正实解不是解集里的第一个时他就取不了，用上边的for取出的也不行，说是syms的

x_positive = xc(xc>0) %%%%%%%%%%%%%%%正值赋值语句

帮我调一下呀，万分感激

Answer 1

[sipek]

有两点修改：

1. 你只需要正实根，所以，对solve函数可以加一个 'Real' 参数，限定只求实根，

2. 在 x 为实根的前提下，用x > 0来才能取出正的实根（如果x为复数，这个判断并不能排除复数）

syms p1 p2 t0 t2 u K L fai2 bai xm zz z
%%%%%%%%%%%%%%%%%%%%%%
K=1.4;u=(K-1)/(K+1);
%%%%%%%%%       P2点                                  %%%%%%%
px2=3;py2=3;pz2=0;Q=1;
%%%%%%%%%%%%%%%%%%%%%%%%   爆心S坐标
sx2=9;sy2=3;sz2=1.75;p0=0.1013;ca=340;
%%%%%%%%%%%%%%%%%%%%%%%%%%%
R2=sqrt((sx2-px2)^2+(sy2-py2)^2+(sz2-pz2)^2)   %%%%%%%%%%%%%%% 入射距离
Z2=R2/Q^(1/3)
%%%%%%%%%%%%%%%%%%%%%   入射角
fai2=acos(sz2/R2)*180/pi
%%%%%%%%%%%%%%%%%%%
p_i2=0.084/Z2+0.27/Z2^2+0.7/Z2^3  
pai1_2=(p_i2+p0)/p0;      
%%%%%%%%%%%%%%%%简化
Mi=sqrt(p_i2/(p0*(1+u))+1)
bai_JX2=1.75/(Mi-1)+39
%%%%%%%%%%%%%%%%简化
%%%%%%%%%%%%%%计算值
M=(pai1_2-1)*t0/(pai1_2+u+(1+pai1_2*u)*t0^2);
arfa=-M^2*(1-u)^2+2*M*(1-u)*t0+1;
beta=-M^2*(1-u)^2*t0-M*(1-2*u-t0^2)+t0;
tuga=M*(1+t0^2);
y=arfa^2-4*beta*tuga;
%%y=eval(x^3-Bc*x^2/Ac-Cc*x/Ac-Dc/Ac) %%赋值语句
x = solve(y,'Real',1)
x_positive = x(x>0)

另外，如果你不希望是一长串符号表达式，可以将最后一句换成

x_positive = double(x(x>0))