matlab 立方插值介绍,三次插值函数的Matlab实现

科学计算课的上机作业，留下来供日后参考。

mSpline文件

clear,clc

​

%定义全局变量，方便在函数ms中调用

​

global x;

​

global y;

​

global M;

​

global n;

​

global df0;

​

global dfn;

​

%输入数据

​

x=[0:10];

​

y=[2.51,3.30,4.04,4.70,5.22,5.54,5.78,5.40,5.57,5.70,5.80];

​

df0=0.8;

​

dfn=0.2;

​

%点的个数

​

nt=size(x);

​

n=nt(2)-1;

​

%计算h,f

​

h=zeros(1);

​

f=zeros(1);

​

for k=1:n

​

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

f(k)=(y(k+1)-y(k))./h(k);

​

end

​

%计算lambda,mu,e

​

lambda=zeros(1);

​

mu=zeros(1);

​

me=zeros(1);

​

for k=2:n

​

lambda(k)=h(k)./(h(k-1)+h(k));

mu(k)=h(k-1)./(h(k-1)+h(k));

me(k)=3*(lambda(k)*f(k-1)+mu(k)*f(k));

​

end

​

%构造矩阵AM=B

​

A=zeros(n-1,n-1);

​

A(1,1)=2;

​

A(1,2)=mu(2);

​

for i=2:n-2

​

A(i,i-1)=lambda(i+1);

A(i,i)=2;

A(i,i+1)=mu(i+1);

​

end

​

A(n-1,n-2)=lambda(n);

​

A(n-1,n-1)=2;

​

M=zeros(n-1,1);

​

B=zeros(n-1,1);

​

B(1)=me(2)-lambda(2)*df0;

​

for i=2:n-2

​

B(i)=me(i+1);

​

end

​

B(n-1)= me(n) - mu(n)*dfn;

​

%求解M

​

M=A\B;

​

mvx=0:0.01:10;

​

%调用ms生成插值函数

​

mvy=0;

​

for i=1:1001

​

mvy(i)=ms(mvx(i));

​

end

​

plot(x,y,'o',mvx,mvy,'black');

​

legend('原始数据','三次样条插值')

ms.m文件

function z = ms(mx)

%MS 计算s

​

%声明全局变量

global M;

global x;

global y;

global n;

global df0;

global dfn;

​

%这是添加了初始条件后的M矩阵

mM=[df0;M;dfn];

​

%这里的k表示书中的k+1

z=0;

for k=1:n+1

z=z+(y(k)*malpha(k,mx)+mM(k)*mbeta(k,mx));

end

​

end

​

%实现的alpha函数

function my=malpha(k,mx)

global x;

global y;

global n;

​

if k==1

if x(k)<= mx & mx <= x(k+1)

my=(1+2*(mx-x(k))/(x(k+1)-x(k)))*((mx-x(k+1))/(x(k)-x(k+1)))^2;

else

my=0;

end

return

end

​

if 2<= k & k<= n

if x(k-1) <= mx & mx <= x(k)

my=(1+2*(mx-x(k))/(x(k-1)-x(k)))*((mx-x(k-1))/(x(k)-x(k-1)))^2;

else

if x(k) < mx & mx <= x(k+1)

my=(1+2*(mx-x(k))/(x(k+1)-x(k)))*((mx-x(k+1))/(x(k)-x(k+1)))^2;

else

my=0;

return

end

end

end

​

if k==n+1

if x(1)<= mx & mx < x(k-1)

my=0;

else

if x(k-1) <= mx & mx <= x(k)

my=(1+2*(mx-x(k))/(x(k-1)-x(k)))*((mx-x(k-1))/(x(k)-x(k-1)))^2;

else

my=0;

return

end

end

end

end

​

function my=mbeta(k,mx)

global x;

global y;

global n;

​

if k==1

if x(k)<= mx & mx <= x(k+1)

my=(mx-x(k))*((mx-x(k+1))/(x(k)-x(k+1)))^2;

else

my=0;

end

return

end

​

if 2<= k & k<= n

if x(k-1) <= mx & mx <= x(k)

my=(mx-x(k))*((mx-x(k-1))/(x(k)-x(k-1)))^2;

else

if x(k) < mx & mx <= x(k+1)

my=(mx-x(k))*((mx-x(k+1))/(x(k)-x(k+1)))^2;

else

my=0;

return

end

end

end

​

if k==n+1

if x(1)<= mx & mx < x(k-1)

my=0;

else

if x(k-1) <= mx & mx <= x(k)

my=(mx-x(k))*((mx-x(k-1))/(x(k)-x(k-1)))^2;

else

my=0;

return

end

end

end

end