matlab中positive,matlab中出现Difference order N mus

该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

syms p q r s t k R

syms u i j a

syms x y w

syms cx cy

syms p1

h=35;

l=11;

x=-30:1:30;

y=-30:1:30;

[X,Y]=meshgrid(x,y);

m=size(X);

Z=ones(m);

P1=ones(m);

P2=ones(m);

j=1;

for x1=-30:1:30

i=1;

for y1=-30:1:30

d=x-h/2+l;

g=(1/2).*(d+((y.^2+h^2/4)./d));

if x1<6.5

w=h/2-l+g+(g.^2-h^2/4).^(1/2);

elseif x1==6.5

w=x;

else

w=h/2-l+g-(g.^2-h^2/4).^(1/2);

end

u=w;

%f=1/84+(1/63-1/84).*((u+11)/35);

%m n取值后

f=1/84+(1/63-1/84).*(70*((u+l)./h).^4-224*((u+l)./h).^5+280*((u+l)./h).^6-160*((u+l)./h).^7+35*((u+l)./h).^8);

%J=(u.*(u + 232))./17640;

%m n取值后

J=(u.*(u.^8 - 81*u.^7 + 1116*u.^6 + 60564*u.^5 - 1318464*u.^4 - 27143424*u.^3 + 562028544*u.^2 + 15897378816*u + 1865020561164))./145921525312500;

theta=asin(J); %theta

p=cos(theta);

%p=simple(p); %cos

n=tan(theta);

%n=simple(n); %tan

%任一点球心坐标

%任一点球心坐标

fx=u-(1./f).*J;

fy=0;

%fz=(1./f).*p+(0.00002951*u.^3+0.006562*u.^2-0.0009396*u+0.0008536);

%Cn=1,拟合曲线取6阶

%fz=(1./f).*p+(0.00000000007519*u.^6+0.000000002433*u.^5+0.0000002793*u.^4+0.00001898*u.^3+0.006576*u.^2-0.000002726*u-0.0000009864);

%m n取值后

fz=(1./f).*p+(2.277e-07*u.^4 + 3.704e-05*u.^3 +0.006494*u.^2 -0.0002022.*u-0.00154);

Z=fz-((1./f).^2-(x-fx).^2-y.^2).^(1/2);

%%%%柱面球面表达

P1=diff(Z,'x');

%q=diff(Z,'y');

%s=diff(p1,'y');

P2=subs(P1,{x,y},{cx,cy});

%q1=subs(q,{x,y},{x1,y1});

%s1=subs(s,{x,y},{x1,y1});

i=i+1;

end

j=j+1;

end