clc;clear;
syms t t1 Eds
%%给定弓高误差,求离散曲线步长
Err=0.01;
tr=0.61166249619097100371595901319577;
hold on
%%绘制分割曲线图案
for tx=0.1:0.01:0.7
x=cos(tx)*sin(tx/2);
y=sin(tx)*sin(tx/2);
plot(x,y,'o');
end
hold off
x=10*cos(t)*sin(t/2);
y=10*sin(t)*sin(t/2);
x1=10*cos(t1)*sin(t1/2);
y1=10*sin(t1)*sin(t1/2);
L=quadl(inline('(5.*2.^(1/2).*(5 - 3.*cos(t)).^(1/2))/2'),0.1,0.7)%%计算弧长
dx=simplify(diff(x,t));
dy=simplify(diff(y,t));
ddx=simplify(diff(dx,t));
ddy=simplify(diff(dy,t));
dx1=simplify(diff(x1,t1));
dy1=simplify(diff(y1,t1));
ddx1=simplify(diff(dx1,t1));
ddy1=simplify(diff(dy1,t1));
dd=[ddx,ddy];
dd1=[ddx1,ddy1];
f=norm(dd,2)%%2范数
f1=norm(dd1,2)%%2范数
ff=-f;
ft=sqrt((8*Err)/f);%%给的一个初值,由于ddc是按初值给的,可能按这个初值去求dt不合适
%%会产生误差,会在该ti+1和ti之间的误差超差,即得在区间内重新求出最大的c的范数对应tmax,
%%这时tmax为该误差下分割最大的tmax,在重新按误差算出dt与ti相加小于tmax即可,不然只能dt0=dt1
ft1=sqrt((8*Err)/f1);
t0=0.1;
t=t0;
dt=eval(ft);
j=1;
while 1
ti=t0+dt;
ta=double(t0);
tb=double(ti);
goalf=inline(ff);
tmax=fminbnd(goalf,ta,tb);%%求函数的最小值,即求曲线的曲率半径最小对应的t值,和其R值
t=tmax;
goalC=eval(f);
goalt=sqrt((8*Err)/goalC);
if (ta+goalt)>tmax
sd(j)=tb;
fd(j)=quadl(inline('(5.*2.^(1/2).*(5 - 3.*cos(t)).^(1/2))/2'),ta,tb);
t0=tb;
t=t0;
ft=sqrt((8*Err)/f);
dt=eval(ft);
else
dt=goalt;
sd(j)=t0+dt;
fd(j)=quadl(inline('(5.*2.^(1/2).*(5 - 3.*cos(t)).^(1/2))/2'),ta,tb);
t0=t0+dt;
ft=sqrt((8*Err)/f);
dt=eval(ft);
end
if abs(sd(max(size(sd)))-0.7)>eps&&j==7
ft1=sqrt((8*Err)/f1);
eqn=0.7==t1+ft1;
tr=solve(eqn,t1)
Err2=sd(max(size(sd))-1)-tr;
Err=(1-Err2)*Err;
t0=0.1;
t=t0;
ft=sqrt((8*Err)/f);
dt=eval(ft);
if Err2<eps
break;
end
end
if sd(max(size(sd)))<0.7&&j==7
break;
end
if j==7
j=0;
sd=[];
fd=[];
end
j=j+1;
end
等精度分割曲线(问题篇)
最新推荐文章于 2021-03-22 01:37:21 发布