MATLAB三次样条插值之三弯矩法

24.构造三次样条函数s(x)去模拟一只飞鸟外形的上部，测得的数据如下：

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,1 1.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 重新输入');

end

v=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);

end

v(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));

end

d(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 重输');

end

h=zeros(n,1);

h(1)=x(1)-x0;

for k=2:n

h(k)=x(k)-x(k-1);

end

[a,b,c,d]=sanzhi(x,y,x0,y0,y1a,y1b);

%追赶解矩阵%

m=chase(a,b,c,d);

M=[1;m;0];

fprintf('三次样条插值的表达式\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);

end

S=S.';

disp(S);

fprintf('习惯解析式如下：\n');

pretty(S);

以下为飞鸟上部数据图形，三次样条曲线图形：

图1 程序运行结果图

图2 利用Matlab的样条差值函数得到的结果图