clc;clear;
syms t Eds
x=10*cos(t)*sin(t/2);
y=10*sin(t)*sin(t/2);
dx=simplify(diff(x,t));
dy=simplify(diff(y,t));
quadl('Fcure(t)',0.1,0.7)
L=quadl(inline('(5.*2.^(1/2).*(5 - 3.*cos(t)).^(1/2))/2'),0.1,0.7)%%计算弧长
r=simplify(sqrt(dx^2+dy^2));
hold on
%%绘制需要分割的图形
for t=0.1:0.01:0.7
x=cos(t)*sin(t/2);
y=sin(t)*sin(t/2);
plot(x,y,'o-');
end
%%按0.2弧长长度进行划分
Eds=0.2;
%%基于0.2划分出来的段数
Ni=L/Eds
%%根据四舍五入取正,那么问题来了,我们知道这样做曲线的最后一段的弧长误差过大,如何解决,采用等弧长划分
%%用整段的实际弧长除以(0.2乘以N(16)得出比例),等同于按Ni除以N,算出dA,即我们提出目标,如何取t
%%可以把产生的误差,通过新的分割取消掉,方法通过牛顿迭代法找出所以的t值,即得出了新的分割长度,当然有人会问
%%都知道比例了,直接0.2*dA不就算出分割新的弧长了嘛,哈哈,好像就是这样求的.
N=16;
Err=N*0.20-L;%%分割前的误差
dA=Ni/N;
t0=eps+0.1;
S=quadl(inline('(5.*2.^(1/2).*(5 - 3.*cos(t)).^(1/2))/2'),0.1,t0);
ft0=S/Eds%%计算弧长
ft0=ft0-dA;%%构造牛顿迭代法
dS=((5*2^(1/2)*(5 - 3*cos(t0))^(1/2))/2);%%构造函数的求导
dft0=dS/Eds
j=0;
ta=0.1;
while 1
ti=t0-ft0/dft0;
fti=quadl(inline('(5.*2.^(1/2).*(5 - 3.*cos(t)).^(1/2))/2'),ta,ti)/Eds%%计算弧长
fti=fti-dA;
dfti=((5*2^(1/2)*(5 - 3*cos(ti))^(1/2))/2)/Eds
if abs(ti-t0)<eps||abs(fti)<eps
j=j+1;
st(j)=ti;
ft(j)=quadl(inline('(5.*2.^(1/2).*(5 - 3.*cos(t)).^(1/2))/2'),0.1,ti);
ta=ti;
t0=ta+eps;
S=quadl(inline('(5.*2.^(1/2).*(5 - 3.*cos(t)).^(1/2))/2'),ta,t0);
fti=S/Eds%%计算弧长
fti=fti-dA;%%构造牛顿迭代法
dS=((5*2^(1/2)*(5 - 3*cos(t0))^(1/2))/2);%%构造函数的求导
dfti=dS/Eds
end
t0=ti;
ft0=fti;
dft0=dfti;
if j==N
break;
end
end
Eds=ft(2)-ft(1)%%新的分割
Err=Eds*N-L%%分割后的误差