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

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=nerror(x or y 输入有误，再来);endv=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-1h(k)=x(k+1)-x(k);endv(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-1v(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);endnewv=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=nerror(x or y 输入有误，再来);end%重新定义hh=zeros(n,1);h(1)=x(1)-x0;for k=2:nh(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:nfprintf(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);endS=S.;disp(S);fprintf(以上为样条函数(三弯矩)解析式，显示为手写如下：n);pretty(S);%第二部分%是否继续运行程序myloop=input(继续运行程序输入“1”，否则输入“0”n);if myloopwhile myloopxi=input(输入需要计算的点的值，并按回车键n);if xix0|xix(n)fprintf(现在开始计算输入点的插值函数值n);elsefprintf(输入数值不在插值范围内，请重新输入n);xi=input(输入需要计算的点的值，并按回车键n);end% 确定输入的数值应该使用哪个解析式newx=x0;x;r,suoy=min(abs(newx-xi);fprintf(输入点的插值函数值为:nn);fprintf(t);if xi=newx(suoy)f=subs(S(suoy-1),X,xi);elsef=subs(S(suoy),X,xi);enddisp(f);myloop=input(继续计算输入“1”，终止计算输入“0”n);endelsereturn;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-1v(i)=c(i)/(b(i)-v(i-1)*a(i-1);u(i)=(d(i)-u(i-1)*a(i-1)/(b(i)-v(