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

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

首先说这个程序并不完善，为了实现通用（1，2，…，n）格式解题，以及为调用追赶法程序，没有针对节点数在三个以下的情况进行分类讨论。希望能有朋友给出更好的方法。

首先，通过函数 sanwanj得到方程的系数矩阵，即追赶法方程的四个向量参数，接下来调用

追赶法（在intersanwj函数中），得到三次样条分段函数系数因子，然后进行多项式合并得

到分段函数的解析式，程序最后部分通过判断输入值的区间自动选择对应的分段函数并计算改

点的值。附：追赶法程序 chase

%%%%%%%%%%%%%%

function [newv,w,newu,newd]=sanwj(x,y,x0,y0,y1a,y1b)

% 三弯矩样条插值

% 将插值点分两次输入，x0 y0 单独输入

% 边值条件a的二阶导数 y1a 和b的二阶导数 y1b，这里建议将y1a和y1b换成y2a和

y2b，以便于和三转角代码相区别

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

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)=6*((y(2)-y(1))/h(1)-(y(1)-y0)/h0)/(h0+h(1));

%

for k=2:n-1

v(k)=h(k-1)/(h(k-1)+h(k));

u(k)=1-v(k);

d(k)=6*((y(k+1)-y(k))/h(k)-(y(k)-y(k-1))/h(k-1))/(h(k-1)+h(k));

end

newv=[v;1];

newu=[1;u];

d0=6*((y(1)-y0)/h0-y1a)/h0;

d(n)=6*(y1b-(y(n)-y(n-1))/h(n-1))/h(n-1);

newd=[d0;d];

%%%%%%%%%%%%

function intersanwj(x,y,x0,y0,y1a,y1b)

% 三弯矩样条插值

%第一部分

n=length(x);m=length(y);

if m~=n

error('x or y 输入有误，再来');

end

%重新定义h

h=zeros(n,1);

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

for k=2:n

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

end

%sptep1 调用三弯矩函数

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

% 三对角方程

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

% 求插值函数

fprintf('三次样条(三弯矩)插值的函数表达式\n');

syms X ;

fprintf('S0--1:\n');

S(1)=collect(((1/6)*M(2)*(X-x0).^3-(1/6)*M(1)*(X-x(1)).^3+(y(1)-

(M(2)*h(1).^2)/6)*(X-x0)-(y0-(M(1)*h(1).^2)/6)*(X-x(1)))/h(1));

for k=2:n

fprintf('S%d--%d:\n',k-1,k);

S(k)=collect(((1/6)*M(k+1)*(X-x(k-1)).^3-(1/6)*M(k)*(X-x(k)).^3+(y(k)-(M(k+1)*h(k).^2)/6)*(X-x(k-1))-(y(k-1)-(M(k)*h(k).^2)/6)*(X-x(k)))/h(k)); end

S=S.';

disp(S);

fprintf('以上为样条函数(三弯矩)解析式，显示为手写如下：\n');

pretty(S);

%第二部分

%是否继续运行程序

myloop=input('继续运行程序输入“1”，否则输入“0”\n');

if myloop

while myloop

xi=input('输入需要计算的点的值，并按回车键\n');

if xi>x0|xi

fprintf('现在开始计算输入点的插值函数值……\n');

else

fprintf('输入数值不在插值范围内，请重新输入\n');

xi=input('输入需要计算的点的值，并按回车键……\n');

end

% 确定输入的数值应该使用哪个解析式

newx=[x0;x];

[r,suoy]=min(abs(newx-xi));

fprintf('输入点的插值函数值为:\n\n');fprintf('\t');

if xi<=newx(suoy)

f=subs(S(suoy-1),X,xi);

else

f=subs(S(suoy),X,xi);

end

disp(f);

myloop=input('继续计算输入“1”，终止计算输入“0”\n');

end

else

return;

end

%%%%%%%%%%%%

function [x]=chase(a,b,c,d)

%追赶法解性方程组 a是下三角b是对角线c是上三角 d是常数项%输入的a b c d 均为列向量

n=length(b);

u=zeros(n,1);

v=zeros(n,1);

x=zeros(n,1);

%追

v(1)=c(1)/b(1);u(1)=d(1)/b(1);

for i=2:n-1

v(i)=c(i)/(b(i)-v(i-1)*a(i-1));

u(i)=(d(i)-u(i-1)*a(i-1))/(b(i)-v(i-1)*a(i-1));

end