MATLAB中有自带的三次样条插值函数
 
插值函数:spline,调用方法： yk =spline(x,y,xk)
 例子：
 
% 读取插值数据
a = load('data.txt');
x = a(:,1);
y = a(:,2);
xx = linspace(max(x),min(x),100);
yy = spline(x,y,xx);    %三次样条插值

% 保存插值结果
fileID = fopen('result.txt','w+');
fprintf(fileID,'%4s %12s\r\n','xx','yy');
fprintf(fileID,'%6.2f     %12.8f\n',xx,yy);
fclose(fileID)

%画图
plot(x,y,'o',xx,yy,'-r')
legend('true','cubicSpline')
saveas(gcf,'out.jpg')
 
data.txt
 
1	93
30	96
60	84
90	84
120	48
150	38
180	51
210	57
240	40
270	45
300	50
330	75
360	80
390	60
420	72
450	67
480	71
510	7
540	74
570	63
600	69
 
计算结果
 
 

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

%x,y为原始离散数据横坐标，纵坐标
%xi为插值间隔，本程序设置为0.1
%调用spline函数即可得到插值结果，再绘图即可
x = 0:2:22;
y = 2*exp(x).*sin(x);
xi = 0:0.1:22;
yi = spline(x,y,xi);
plot(x,y,'*',xi,yi);
 
参考资料：https://jingyan.baidu.com/article/fb48e8beda889e6e622e14e6.html

 
 
一些语句
 
 
size(a,x);
 
 
 
  
xis=size(xi,2);  %获取矩阵xi的列数
xis=size(xi,1);  %获取矩阵xi的行数
 
 
 
 
round函数 四舍五入，在绘制三维图形时如果数据比较大，可以通过以下方法把数据按应有比例减小  （我自己理解的，有错请指正）
 
 
 
  
x=round(x'./100+1);  %这部分使数据减小，便于作图处理等
y=round(y'./100+1);
 
 
 
 
matlab中提供的插值函数 griddata ，其基本格式为 cz=griddata(x,y,z,cx,cy,'method');
 
 
其中，参数method有以下四个可以选择：
 
 
‘nearest’ ：最邻近点插值
 
 
‘linear’ ：线性插值
 
 
‘cubic’： 三次插值
 
 
‘V4’：matlab中所提供的的插值方法
 
 
基本格式中 x,y,z均是n为向量，代表数据点的坐标，向量cx cy是给定的网格点的横坐标和纵坐标，cz是矩阵，其行数和列数分别等于 cx 和 cy 的维数。（所以前面有提到round函数的用处，不然cz矩阵可能会特别特别大）
 
 
 
  
%手写一下代码
clear
clc

%读入数据
load textx.txt
load texty.txt
load textz.txt

%数据进行处理
textx=round(textx'./100+1);
texty=round(texty'./100+1);
textz=textz';

zi=zeros(max_x,max_y); %这里max_x max_y我理解的是textx texty中最大的值（初学者）
for i=1:len
     zi(x(i),y(i))=z(i);
end

xi=0:step:lx; %step为步长  lx为自己设定的x轴长度，要大于等于max_x ，ly同理
yi=0:step:ly
zi=griddata(x,y,z,xi,yi,'V4');%绘制三维曲面图
subplot(1,2,1);
mesh(xi,yi,zi); %绘制地形图

%这是最基本的结构，其他内容可以在此基础上添加
 
 
 
 
 
 
 
 
 
 
输出
 
 
转载于:https://www.cnblogs.com/yz-lucky77/p/11305119.html