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

非常类似前面的三弯矩法，这里的sanzhj函数和intersanzhj作用相当于前面的sanwanj和intersanwj，追赶法程序通用，代码如下。%%%%%%%%%%%%%%%%%%%function[newu,w,newv,d]=sanzhj(x,y,x0,y0,y1a,y1b)%三转角样条插值%将插值点分两次输入，x0y0单独输入%边值条件a的一阶导数y1a和b的一阶导数y1bn=length(x);m=length(y);ifm~=nerror('xory输入有误，再来');endv=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);fork=1:n-1h(k)=x(k+1)-x(k);endv(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);%fork=2:n-1v(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));endd(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,:);%%%%%%%%%%%%functionintersanzhj(x,y,x0,y0,y1a,y1b)%三转角样条插值%第一部分n=length(x);m=length(y);ifm~=nerror('xory输入有误，再来');end%重新定义hh=zeros(n,1);h(1)=x(1)-x0;fork=2:nh(k)=x(k)-x(k-1);end%调用三转角函数[a,b,c,d]=sanzhj(x,y,x0,y0,y1a,y1b);%三对角方程m=chase(a,b,c,d);%求MM=[1;m;0];%求插值函数fprintf('三次样条（三转角）插值的函数表达式\n');symsX;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);fork=2:nfprintf('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);endS=S.';disp(S);fprintf('以上为样条函数（三转角）解析式，显示为手写习惯如下：\n');pretty(S);%第二部分%是否继续运行程序myloop=input('继续运行程序输入“1”，否则输入“0”\n');ifmyloopwhilemyloopxi=input('输入需要计算的点的值，并按回车键\n');ifxix0|xix(n)fprintf('现在开始计算输入点的插值函数值……\n');elsefprintf('输入数值不在插值范围内，请重新输入\n');xi=input('输入需要计算的点的值，并按回车键……\n');end%确定输入的数值应该使用哪个解析式newx=[x0;x];[r,suoy]=min(abs(newx-xi));fprintf('输入点的插值函数值为:\n\n');fprintf('\t');ifxi=newx(suoy)f=subs(S(suoy-1),X,xi);elsef=subs(S(suoy),X,xi);enddisp(f);myloop=input('继续计算输入“1”，终止计算输入“0”\n');endelsereturn;end