程序代码如下
function [length]=contactlength(x0)
if x0>50||x0<-50
    'error:数据超出尺寸范围'
else
    syms x
    R=300;%非球面顶点曲率半径
    c=1/R;
    delta=0.1;
    k=-3.3;%非球面参数
    rb=27;%半径
    y=(-1*c*x.^2)./(1+sqrt(1-(1+k)*(c^2)*x.^2));
    dy=diff(y);
    dy2=diff(y,2);
    dyx0=subs(dy,x0);
    dy2x0=subs(dy2,x0);
    r0=(1+dyx0^2)^(3/2)/(abs(dy2x0))%曲率半径
    y0=(-1*c*x0^2)./(1+sqrt(1-(1+k)*(c^2)*x0^2));%非球面接触点y坐标
   if x0==0
      xb=x0;
      yb=y0+rb-delta;
      ob=[xb,yb]%圆心二维坐标
   else
      xb=x0+(rb-delta)/sqrt(1+1/dyx0^2);
      yb=y0+(-1/dyx0)*(rb-delta)/sqrt(1+1/dyx0^2);
      ob=[xb,yb]%圆心二维坐标
   end
    syms x2 y2
    ezplot((x2-xb)^2+(y2-yb)^2==rb^2,[-50 50 -10 70])
    axis equal;
    hold on;
    ezplot(y2==(-1*c*x2^2)/(1+sqrt(1-(1+k)*(c^2)*x2^2)),[-50 50 -10 70])
    s=solve(y2==(-1*c*x2^2)./(1+sqrt(1-(1+k)*(c^2)*x2^2)),...
            (x2-xb)^2+(y2-yb)^2==rb^2,x2,y2);
    X=double(s.x2);
    Y=double(s.y2);
    X=X(~logical(imag(X)));
    length=max(X)-min(X)
end
在实际运算过程中，要100多秒才能算出结果，请问是哪里有问题呢？我该怎么改程序？
本人新手，接触matlab不久，希望大神们多多指导 返回小木虫查看更多