MATLAB三次樣條插值之三彎矩法.doc

24.構(gòu)造三次樣條函數(shù)s(x)去模擬一只飛鳥外形的上部，測得的數(shù)據(jù)如下：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.2,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; 邊界條件為自然邊界條件，即：S(0.9)=0,s(13.3)=0（可畫出圖形）。解：程序：%解樣條差值%function newu,w,newv,d=sanzhi(x,y,x0,y0,y1a,y1b)n=length(x);m=length(y);if m=n error(x or y重新輸入);endv=ones(n-1,1);u=ones(n-1,1);d=zeros(n-1,1);w=2*ones(n-1,1);h0=x(1)-x0;h=zeros(n-1,1);for k=1:n-1 h(k)=x(k+1)-x(k);endv(1)=h0/(h0+h(1);u(1)=1-v(1);d(1)=3*(v(1)*(y(2)-y(1)/h(1)+u(1)*(y(1)-y0)/h0);for k=2:n-1 v(k)=h(k-1)/(h(k-1)+h(k); u(k)=1-v(k); d(k)=3*(v(k)*(y(k+1)-y(k)/h(k)+u(k)*(y(k)-y(k-1)/h(k-1);endd(1)=d(1)-u(1)*y1a;d(n-1)=d(n-1)-v(n-1)*y1b;newv=v(1:n-2,:);newu=u(2:n-1,:);function intersanzhi(x,y,x0,y0,y1a,y1b)n=length(x);m=length(y);if m=n error(x or y重輸);endh=zeros(n,1);h(1)=x(1)-x0;for k=2:n h(k)=x(k)-x(k-1);enda,b,c,d=sanzhi(x,y,x0,y0,y1a,y1b);%追趕解矩陣%m=chase(a,b,c,d);M=1;m;0;fprintf(三次樣條插值的表達(dá)式n);syms X ;fprintf(S0-1:n);S(1)=collect(y0/h(1).3)*(X-x(1).2*(h(1)+2*(X-x0)+(y(1)/h(1).3)*(X-x0).2*(h(1)+2*(x(1)-X)+(M(1)/h(1).2)*(X-x0)*(X-x(1).2+(M(2)/h(1).2)*(X-x(1)*(X-x0).2);for k=2:n fprintf(S%d-%d:n,k-1,k); S(k)=collect(y(k-1)/h(k).3)*(X-x(k).2*(h(k)+2*(X-x(k-1)+(y(k)/h(k).3)*(X-x(k-1).2*(h(k)+2*(x(k)-X)+(M(k)/h(k).2)*(X-x(k-1)*(X-x(k).2+(M(k+1)/h(k).2)*(X-x(k)*(X-x(k-1).2);endS=S.;disp(S);fprintf(習(xí)慣解