弧长离散曲线

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%%分割后的误差