程序代碼如下

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不久，希望大神們多多指導