非常类似前面的三弯矩法,这里的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<x(n)
            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

 

posted on 2008-12-16 08:46  秋波渡  阅读(3981)  评论(0编辑  收藏  举报