MATLAB三次样条插值之三转角法

非常类似前面的三弯矩法，这里的sanzhj函数和intersanzhj作用相当于前面的sanwanj和intersanwj，追赶法程序通用，代码如下。

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

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

%三转角样条插值

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

%边值条件a的一阶导数y1a 和b的一阶导数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 intersanzhj(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

%调用三转角函数

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

%三对角方程

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

%求M

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

%第二部分

%是否继续运行程序

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