三次样条插值自然边界条件

例：已知一组数据点心■….门，编写一程序求解三次样条插值函数勺,满足月（第£）=Hi•>4 1*2**..、门并针对下面一组具体实验数据J?;0.250.30.390.450.53Vi0.50000.54770.62450.67080.7280求解，其中边界条件为、「1.«」「：；）一U.1）三次样条插值自然边界条件源程序：functions=spline3(x,y,dy1,dyn)%又为节点，y为节点函数值，dyl,dyn分别为x=0.25,0.53处的二阶导m=length(x);n=length(y);ifm~=nerror'xor尸输入有误')returnendh=zeros(1,n-1);h(n-1)=x(n)-x(n-1);fork=1:n-2h(k)=x(k+1)-x(k);v(k)=h(k+1)/(h(k+1)+h(k));u(k)=1-v(k);endg(1)=3*(y(2)-y(1))/h(1)-h(1)/2*dy1;g(n)=3*(y(n)-y(n-1))/h(n-1)+h(n-1)/2*dyn;fori=2:n-1g(i)=3*(u(i-1)*(y(i+1)-y(i))/h(i)+v(i-1)*(y(i)-y(i-1))/h(i-1));endfori=2:n-1;A(i,i-1)=v(i-1);A(i,i+1)=u(i-1);endA(n,n-1)=1;A(1,2)=1;A=A+2*eye(n);M=zhuigf(A,g); %调用函数，追赶法求Mfprintf('三次样条(三对角)插值的函数表达式\n');symsX;fork=1:n-1fprintf'S%d--%d:\n',k,k+1);s(k)=(h(k)+2*(X-x(k)))./h(k).八3.*(X-x(k+1)).八2.*y(k)...+(h(k)-2*(X-x(k+1)))./h(k).八3.*(X-x(k)).八2.*y(k+.)..+(X-x(k)).*(X-x(k+1)).八2./h(k)"2*M(k)+(X-x(k+1)).*..(X-x(k))"2./h(k).八2*M(k+1);ends=s.';s=vpa(s,4);%画三次样条插值函数图像fori=1:n-1X=x(i):0.01:x(i+1);st=(h(i)+2*(X-x(i)))./(h(i)八3).*(X-x(i+1))"2.*y(i)...+(h(i)-2.*(X-x(i+1)))./(h(i)八3).*(X-x(i)).八2.*y(i+1)...+(X-x(i)).*(X-x(i+1)).八2./h(i)八2*M(i)+(X-x(i+1)).*...(X-x(i))"2./h(i)八2*M(i+1);plot(x,y,o',X,st);holdonEndplot(x,y);gridon%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%调用的函数：%追赶法functionM=zhuigf(A,g)n=length(A);L=eye(n);U=zeros(n);fori=1:n-1U(i,i+1)=A(i,i+1);endU(1,1)=A(1,1);fori=2:nL(i,i-1)=A(i,i-1)/U(i-1,i-1);U(i,i)=A(i,i)-L(i,i-1)*A(i-1,i);endY(1)=g(1);fori=2:nY(i)=g(i)-L(i,i-1)*Y(i-1);endM(n)=Y(n)/U(n,n);fori=n-1:-1:1M(i)=(Y(i)-A(i,i+1)*M(i+1))/U(i,i);end2)在命令窗口输入*,y,dy1,dyn，得到三次样条函数：>>x=[0.25,0.3,0.39,0.45,0.53];>>y=[0.5,0.5477,0.6245,0.6708,0.7280];>>dy1=0;>>dyn=0;>>s=spline3(x,y,dy1,dyn)运行结果：三次样条(三对角)插值的函数表达式TOC\o"1-5"\h\z-2:-3:-4:-5:.5000*(-3600.+.1600e5*X)*(X-.3000)A2+.5477*(5200.-.1600e5*X)*(X-.2500)A2+394.8*(X-.2500)*(X-.3000)A2+355.2*(X-.3000)*(X-.2500)A2.5477*(-699.6+2743.*X)*(X-.3900)A2+.6245*(1193.-2743.*X)*(X-.3000)A2+109.6*(X-.3000)*(X-.3900)A2+96.77*(X-.3900)*(X-.3000)A2.6245*(-3333.+9259.*X)*(X-.4500)A2+.6708*(4444.-9259.*X)*(X-.3900)A2+217.7*(X-.3900)*(X-.4500)A2+207.6*(X-.4500)*(X-.3900)A2.6708*(-1602.+3906.*X)*(X-.5300)A2+.7280*(2227.-3906.*X)*(X-.4500)A2+116.8*(X-.4500)*(X-.5300)A2+109.2*(X-.5300)*(X-.4500)A2如将三次样条函数加以整理，可用如下程序：s=collect(s);则输出结果为：s=.4595000000000-13.200*X八3+9.9000000*X八2-1.48800000000*X.37459306800000-4.2924*3+3.66332200*2-.137