三次样条函数模拟飞鸟外形上部(共3页).doc
精选优质文档-倾情为你奉上clear clc x=0.9 1.3 1.9 2.1 2.6 3.0 3.9 4.4 4.7 5.0 6.0 7.0 8.0 9.2 10.5 11.3 11.6 12.0 12.6 13.0 13.3; y=1.3 1.5 1.85 2.1 2.6 2.7 2.4 2.15 2.05 2.1 2.25 2.3 2.25 1.95 1.4 0.9 0.7 0.6 0.5 0.4 0.25; n=length(x); for i=1:n-1 h(i)=x(i+1)-x(i); endfor i=1:n-2 k(i)=h(i+1)/(h(i)+h(i+1); u(i)=h(i)/(h(i)+h(i+1); endfor i=1:n-2 gl(i)=3*(u(i)*(y(i+2)-y(i+1)/h(i+1)+k(i)*(y(i+1)-y(i)/h(i);endg0=3*(y(2)-y(1)/h(1); g00=3*(y(n)-y(n-1)/h(n-1); g=g0 gl g00;g=transpose(g); k1=k 1; u1=1 u; Q=2*eye(21)+diag(u1,1)+diag(k1,-1); m=transpose(Qg);syms X;for i=1:n-1 p1(i)=(1+2*(X-x(i)/h(i)*(X-x(i+1)/h(i)2*y(i); p2(i)=(1-2*(X-x(i+1)/h(i)*(X-x(i)/h(i)2*y(i+1); p3(i)=(X-x(i)*(X-x(i+1)/h(i)2*m(i); p4(i)=(X-x(i+1)*(X-x(i)/h(i)2*m(i+1); p(i)=p1(i)+p2(i)+p3(i)+p4(i); p(i)=simplify(p(i); ends1=p(1) s2=p(2) s3=p(3) s4=p(4) for k=1:20 for z=x(k):0.001:x(k+1) q=eval(subs(p(k),'X','z'); plot(z,q,'b') hold on endendgrid on title('the curve of a bird') xlabel('x') ylabel('p')结果s1 = - (*X3)/ + (*X2)/ - (*X)/ + / s2 = (*X3)/ - (*X2)/ + (*X)/ - / s3 = - (*X3)/70496 + (*X2)/40992 - (*X)/ + / s4 = - (40237*X3)/85248 + (*X2)/13120 - (*X)/ + / 专心-专注-专业
