程序代码如下

function [length]=contactlength(x0)

if x0>50||x0

'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不久，希望大神们多多指导