我已经把东西放在附件里,高手可以放到自己机子上运行一下.谢谢.


R1=0.5*(str2double(get(handles.R1_input,'string')));
R2=0.5*(str2double(get(handles.R2_input,'string')));
R3=0.5*(str2double(get(handles.R3_input,'string')));
e=str2double(get(handles.e_input,'string'));
d=str2double(get(handles.d_input,'string'));
w1=str2double(get(handles.w1_input,'string'));
w2=str2double(get(handles.w2_input,'string'));
Xv=str2double(get(handles.Xv_input,'string'));
Xs=Xv;
syms t;
syms n n1 n2 n3;
n1=0;
n2=0;
%工具间隔l1,沿半径方向分成1000份
l1=(R1-R3)/1000;
%算切入的最大直径
d2m=0.3*d;
%定义MMR
MMR=zeros(1,R2);
for R20=1:1:R2
    R21=R20;
    V=0;
    for m=R3:l1:R1
        s1=e^2+R21^2-2*e*R21*cos(w1*t)-m^2;
        t1=solve(s1,'t');
        x1=e*cos(w1*t1)+R21*cos((w2-w1)*t1);
        y1=(-1)*e*sin(w1*t1)+R21*sin((w2-w1)*t1);
         
        s2=e^2+R21^2-2*e*R21*cos(w1*t)-(m+l1)^2;
        t2=solve(s2,'t');
        x2=e*cos(w1*t2)+R21*cos((w2-w1)*t2);
        y2=(-1)*e*sin(w1*t2)+R21*sin((w2-w1)*t2);
        
        pq=sqrt((x1-x2)^2+(y1-y2)^2);
        %算面积s
        s=pq*2*pi*0.001*Xs;
        %算磨粒数目n;
        %s3=(d^2)*(n+1)*(2*n+1)-(24*n*s);%%%%??????
        n=solve('(d^2)*(n+1)*(2*n+1)-(24*n*s)','n');
        n11=n(1);
        n22=n(2);
        if isreal(n11)==1 &&n11>100
            n1=n11*0.3;
        elseif isreal(n22)==1 &&n22>100
            n2=n22*0.3;
        else   
            continue
        end
        
        if (n1-n2)>0
            n3=n1;
        else
            n3=n2;
        end
        
        %定义切削累加体积v
        v=0;
        for f=1:1:n3
              d2=d2m/n3*f;
              %算切入面积(A1-A2)
              %角度u
              u=asin(d2/d);
              A1=(d/2)^2*u;
              A2=(d2/2)^2*(sin(2*u))*0.5;
              v1=(A1-A2)*d2;
              v=v1+v;
        end
    V=v+V;
    end
    V0=V*R21*1000/R2;
    MMR(R21)=V0;
end
axes(handles.axes1)
plot(MMR,R20) ;        


为什么老有错误??[em2][em2][em2]