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

《MATLAB三次样条插值之三弯矩法》由会员上传分享，免费在线阅读，更多相关内容在应用文档-天天文库。

1、MATLAB三次样条插值之三弯矩法首先说这个程序并不完善，为了实现通用（1，2，…，n）格式解题，以及为调用追赶法程序，没有针对节点数在三个以下的情况进行分类讨论。希望能有朋友给出更好的方法。首先，通过函数sanwanj得到方程的系数矩阵，即追赶法方程的四个向量参数，接下来调用追赶法（在intersanwj函数中），得到三次样条分段函数系数因子，然后进行多项式合并得到分段函数的解析式，程序最后部分通过判断输入值的区间自动选择对应的分段函数并计算改点的值。附：追赶法程序chase%%%%%%%%%%%%%%function[newv,w,newu,newd]=s

2、anwj(x,y,x0,y0,y1a,y1b)%三弯矩样条插值%将插值点分两次输入，x0y0单独输入%边值条件a的二阶导数y1a和b的二阶导数y1b，这里建议将y1a和y1b换成y2a和y2b，以便于和三转角代码相区别n=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);h0=x(1)-x0;h=zeros(n-1,1);fork=1:n-1h(k)=x(k+1)-x(k);endv(1)=h0

3、/(h0+h(1));u(1)=1-v(1);d(1)=6*((y(2)-y(1))/h(1)-(y(1)-y0)/h0)/(h0+h(1));%fork=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=[

4、d0;d];%%%%%%%%%%%%functionintersanwj(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%sptep1调用三弯矩函数[a,b,c,d]=sanwj(x,y,x0,y0,y1a,y1b);%三对角方程M=chase(a,b,c,d);%求插值函数fprintf('三次样条(三弯矩)插值的函数

5、表达式');symsX;fprintf('S0--1:');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));fork=2:nfprintf('S%d--%d:',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、/6)*(X-x(k-1))-(y(k-1)-(M(k)*h(k).^2)/6)*(X-x(k)))/h(k));endS=S.';disp(S);fprintf('以上为样条函数(三弯矩)解析式，显示为手写如下：');pretty(S);%第二部分%是否继续运行程序myloop=input('继续运行程序输入“1”，否则输入“0”');ifmyloopwhilemyloopxi=input('输入需要计算的点的值，并按回车键');ifxi>x0

7、xi

8、输入数值不在插值范围内，请重新输入');xi=input('输入需要计算的点的值，并按回车键……');end%确定输入的数值应该使用哪个解析式newx=[x0;x];[r,suoy]=min(abs(newx-xi));fprintf('输入点的插值函数值为:');fprintf('t');ifxi<=newx(suoy)f=subs(S(suoy-1),X,xi);elsef=subs(S(suoy),X,xi);enddisp(f);myloop=input('继续计算输入“1”，终止计算输入“0”');endelsereturn;e

9、nd%%%%%%%%%%%%funct