科学计算课的上机作业，留下来供日后参考。
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
