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

1、例：已知一组数据点，编写一程序求解三次样条插值函数满足 并针对下面一组具体实验数据0.250.30.390.450.530.50000.54770.62450.67080.7280求解，其中边界条件为.1）三次样条插值自然边界条件源程序：function s=spline3(x,y,dy1,dyn)%x为节点，y为节点函数值，dy1，dyn分别为x=0.25,0.53处的二阶导m=length(x);n=length(y);if m=n error(x or y输入有误) returnendh=zeros(1,n-1);h(n-1)=x(n)-x(n-1);for k=1:n-2 h(k)=x

2、(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;for i=2:n-1 g(i)=3*(u(i-1)*(y(i+1)-y(i)/h(i)+v(i-1)*(y(i)-y(i-1)/h(i-1);endfor i=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); %

3、调用函数，追赶法求Mfprintf(三次样条（三对角）插值的函数表达式n);syms X;for k=1:n-1 fprintf(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+1). +(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); %画三次样条插值函数图像for i=1:n-1 X=x(i

4、):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); hold on End plot(x,y); grid on%调用的函数：%追赶法function M=zhuigf(A,g)n=length(A);L=eye(n);U=zeros(n);for i=1

5、:n-1 U(i,i+1)=A(i,i+1);endU(1,1)=A(1,1);for i=2:n L(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);for i=2:n Y(i)=g(i)-L(i,i-1)*Y(i-1);endM(n)=Y(n)/U(n,n);for i=n-1:-1:1 M(i)=(Y(i)-A(i,i+1)*M(i+1)/U(i,i);end2)在命令窗口输入x，y,dy1,dyn,得到三次样条函数：x=0.25,0.3,0.39,0.45,0.53;y=0.5,0.547

6、7,0.6245,0.6708,0.7280;dy1=0;dyn=0;s=spline3(x,y,dy1,dyn)运行结果：三次样条（三对角）插值的函数表达式S1-2:S2-3:S3-4:S4-5: .5000*(-3600.+.1600e5*X)*(X-.3000)2+.5477*(5200.-.1600e5*X)*(X-.2500)2+394.8*(X-.2500)*(X-.3000)2+355.2*(X-.3000)*(X-.2500)2 .5477*(-699.6+2743.*X)*(X-.3900)2+.6245*(1193.-2743.*X)*(X-.3000)2+109.6*(X

7、-.3000)*(X-.3900)2+96.77*(X-.3900)*(X-.3000)2 .6245*(-3333.+9259.*X)*(X-.4500)2+.6708*(4444.-9259.*X)*(X-.3900)2+217.7*(X-.3900)*(X-.4500)2+207.6*(X-.4500)*(X-.3900)2 .6708*(-1602.+3906.*X)*(X-.5300)2+.7280*(2227.-3906.*X)*(X-.4500)2+116.8*(X-.4500)*(X-.5300)2+109.2*(X-.5300)*(X-.4500)2如将三次样条函数加以整理，可用如下程序：s=collect(s);则输出结