MATLAB.zip_MATLAB求控制点_求曲线误差_翼型离散点反求bezier曲线控制点_贝塞尔 拟合_贝塞尔拟合

num1=xlsread('C:\Users\13475\Documents\MATLAB\naca0012(1).xlsx'); num2=xlsread('C:\Users\13475\Documents\MATLAB\naca0012(2).xlsx'); ex1=num1(:,1); ex2=num2(:,1); ey1=num1(:,2); ey2=num2(:,2); ex1=ex1'; ex2=ex2'; ey1=ey1'; ey2=ey2'; s1=size(ex1,2); s2=size(ex2,2); deltat1=1/(s1-1); t1=0:deltat1:1; a1=ones(s1,1); deltat2=1/(s2-1); t2=0:deltat2:1; a2=ones(s2,1); n=3;% 可以改变的阶数 A1=zeros(n+1); A1(1,1)=1; A1(n+1,n+1)=1; for i=2:n for u=1:n+1 A1(i,u)=sum(nchoosek(n,u-1)*nchoosek(n,i-1)*(a1'-t1).^(n-i+1).*t1.^(i-1).*(a1'-t1).^(n-u+1).*t1.^(u-1)); end end A2=zeros(n+1); A2(1,1)=1; A2(n+1,n+1)=1; for i=2:n for u=1:n+1 A2(i,u)=sum(nchoosek(n,u-1)*nchoosek(n,i-1)*(a2'-t2).^(n-i+1).*t2.^(i-1).*(a2'-t2).^(n-u+1).*t2.^(u-1)); end end Cx1=zeros(n+1,1); Cx1(1,1)=ex1(1,1); Cx1(n+1,1)=ex1(1,s1); for i=2:n Cx1(i)=sum(nchoosek(n,i-1)*ex1.*(a1'-t1).^(n-i+1).*t1.^(i-1)); end Cx2=zeros(n+1,1); Cx2(1,1)=ex2(1,1); Cx2(n+1,1)=ex2(1,s2); for i=2:n Cx2(i)=sum(nchoosek(n,i-1)*ex2.*(a2'-t2).^(n-i+1).*t2.^(i-1)); end Cy1=zeros(n+1,1); Cy1(1)=ey1(1); Cy1(n+1)=ey1(s1); for i=2:n Cy1(i)=sum(nchoosek(n,i-1)*ey1.*(a1'-t1).^(n-i+1).*t1.^(i-1)); end Cy2=zeros(n+1,1); Cy2(1)=ey2(1); Cy2(n+1)=ey2(s2); for i=2:n Cy2(i)=sum(nchoosek(n,i-1)*ey2.*(a2'-t2).^(n-i+1).*t2.^(i-1)); end a1=0.1; U1=inv(A1)*Cx1; V1=inv(A1)*Cy1; U1=[0;U1]; V1=[V1(1);a1*V1(2);V1(2:n+1)]; a2=-0.5; U2=inv(A2)*Cx2; V2=inv(A2)*Cy2; U2=[0;U2]; V2=[V2(1);a2*V2(2);V2(2:n+1)]; %plot(U,V,'r.'); syms t; n=4; x(t)=vpa((1-t)^(n)*t*U1(1),3); y(t)=vpa(t*(1-t)^(n)*V1(1),3); for i=2:n+1 x(t)=x+nchoosek(n,i-1)*(1-t)^(n-i+1)*t^(i-1)*U1(i); x=vpa(x,3); y(t)=y+nchoosek(n,i-1)*(1-t)^(n-i+1)*t^(i-1)*V1(i); y=vpa(y,3); end h(t)=(1-t)^(n)*t*U2(1); k(t)=t*(1-t)^(n)*V2(1); for i=2:n+1 h(t)=h+nchoosek(n,i-1)*(1-t)^(n-i+1)*t^(i-1)*U2(i); h=vpa(h,3); k(t)=k+nchoosek(n,i-1)*(1-t)^(n-i+1)*t^(i-1)*V2(i); k=vpa(k,3); end q1=0; for i=2:s1-1 j=vpa(fzero(@(t) x(t)-ex1(i),[0,1]),3); q1=[q1;vpa(abs(vpa(y(j),3)-ey1(i))/ey1(i),3)]; %q1=vpa(q1,3); end q1; ezplot(x,y,[0,1]); hold on ezplot(h,k,[0,1]); plot(ex1,ey1,'r'); plot(U1,V1); plot(ex2,ey2,'r'); plot(U2,V2); hold off