matlab中调用追赶法求解系数矩阵的三次样条插值

已知数表x1245y1342边界条件S''(x0)=0，S''(x3)=0插值点x=3---------------------------------------------------------分--割--线---------------------------------------------------------使用三转角方程边界条件为第二类边界条件---------------------------------------------------------分--割--线---------------------------------------------------------先写好追赶法的程序把followup.m存入工作路径functionx=followup(a,b,c,d)n=length(d);a(1)=0;%“追”的过程L(1)=b(1);y(1)=d(1)/L(1);u(1)=c(1)/L(1);fori=2:(n-1)L(i)=b(i)-a(i)*u(i-1);y(i)=(d(i)-y(i-1)*a(i))/L(i);u(i)=c(i)/L(i);endL(n)=b(n)-a(n)*u(n-1);y(n)=(d(n)-y(n-1)*a(n))/L(n);%“赶”的过程x(n)=y(n);fori=(n-1):-1:1x(i)=y(i)-u(i)*x(i+1);end---------------------------------------------------------分--割--线---------------------------------------------------------用matlab编程如下:function[s,y0]=spline3(x,y,x0)%x,y为数表x0为插值点s表示插值函数y0为x0对应的插值函数值symstn=length(x);%得出nfori=1:n-1;h(i)=x(i+1)-x(i);endfori=2:n-1;lamda(i)=h(i)/(h(i-1)+h(i));miu(i)=1-lamda(i);g(i)=3*(lamda(i)*((y(i)-y(i-1))/h(i-1))+miu(i)*((y(i+1)-y(i))/h(i)));endg(1)=3*((y(2)-y(1))/h(1));g(n)=3*((y(n)-y(n-1))/h(n-1));%前边求出lamdamiu和g从而可以确定系数矩阵miu(1)=1;miu(4)=0;lamda(n)=1;%根据第二边界条件补充两个lamda和miu的值fori=1:nbeta(i)=2;endm=followup(lamda,beta,miu,g)%解出m的值从而可确定stst为各段的插值多项式fori=1:n-1st(i)=(t-x(i+1))^2*(h(i)+2*(t-x(i)))*y(i)/(h(i)^3)...+(t-x(i))^2*(h(i)+2*(x(i+1)-t))*y(i+1)/(h(i)^3)...+(t-x(i))^2*(t-x(i+1))*m(i+1)/(h(i)^2)...+(t-x(i+1))^2*(t-x(i))*m(i)/(h(i)^2);end%得到插值的结果各段的t的表达式%接下来要将插值点x0代入首先确定x0所在的插值区间fori=1:n-1if(x(i)x0)in=i;endends=st(in);s=expand(s);s=collect(s,'t');y0=subs(s,'t',x0)%s是插值多项式y0是插值点的函数值---------------------------------------------------------分--割--线---------------------------------------------------------在matlab中输入x=[1245]；y=[1342]；spline3(x,y,3)---------------------------------------------------------分--割--线---------------------------------------------------------得到y0=4.2500ans=-1/8*t^3+3/8*t^2+7/4*t-1