紧压三次样条插值程序讲解.docVIP

求压紧三次样条函数的函数程序 function S liti05_7 X,Y,dx0,dxn %Input -X is the 1xn abscissa vector % -Y is the 1xn ordinate vector % -dx0 S x0 first derivative boundary condition % -dxn S xn first derivative boundary condition %Output -S:rows of S are the confficients,indescending order,for the cubic interplants N length X -1; %求当前问题的规模数--的最大下标 H diff X ; %求X的差分 h0 h1 …… hn-1 D diff Y ./H; %求 y 对 x 的一阶差商 d0 d1 …… dn-1 A H 2:N-1 ; %为求线性方程组系数矩阵上次对角元做准备 B 2* H 1:N-1 +H 2:N ; %求线性方程组系数矩阵主对角元做准备 C H 2:N-1 ; %求线性方程组系数矩阵上次对角元 U 6*diff D ; %为求线性方程组常数列向量做准备 %压紧样条端点约束 B 1 B 1 -H 1 /2; %3h0/2+2h1 U 1 U 1 -3* D 1 -dx0 ; %修改 u 1 B N-1 B N-1 -H N /2; % 2h N-2 +3*h N-1 /2 U N-1 U N-1 -3* dxn-D N ; %修改 u N-1 A %输出线性方程组系数矩阵 上次对角元向量 B %输出线性方程组系数矩阵主对角元向量 C %输出线性方程组系数矩阵下次对角元向量 U %输出线性方程组常数列向量 % 下面开始解三对角线行方程组 % 首先转化为上三角线性方程组 for k 2:N-1 temp A k-1 /B k-1 ; B k B k -temp*C k-1 ; U k U k -temp*U k-1 ; end %回代求解 M N U N-1 /B N-1 ; for k N-2:-1:1 M k+1 U k -C k *M k+2 /B k ; end %计算端点x0 xn上的二阶导数 M 1 3* D 1 -dx0 /H 1 -M 2 /2; M N+1 3* dxn-D N /H N -M N /2; M % 输出各节点上的二阶导数 for k 0:N-1 %计算第k个多项式的系数 S k+1,1 M k+2 -M k+1 / 6*H k+1 ; %算（x-X k ）三次项系数 S k+1,2 M k+1 /2; %算（x-X k ）二次项系数 S k+1,3 D k+1 -H k+1 * 2*M k+1 +M k+2 /6;%算（x-X k ）一次项系数 S k+1,4 Y k+1 ; %算（x- X k ）零次项系数 end hold on plot X,Y,kO % 画样点 style [r, g, b]; % 由于例题中N等于3，style只构造了三个元素 for k 1:N %画第k个区间上的三次函数曲线段 x1 X k :0.01:X k+1 ; y1 polyval S k,: ,x1-X k ; %算关于（x- X k ）三次多项式的值 plot x1,y1,style k end grid on xlabel x ; ylabel y ; hold off 例 求压紧三次样条曲线，经过点（0，0），（1，0.5），（2，2.0），（3，1.5），一阶导数的边界条件为 S’ 0 0.2和S’ 3 -1 针对上例，程序运行结果如下 ?x [0 1 2 3];y [0 0.5 2.0 1.5]; ?dx0 0.2;dxn -1; ?s liti05_7 x,y,dx0,dxn A 1 B 3.5000 3.5000 C 1 U 5.1000 -10.5000 M -0.3600 2.5200 -3.7200 0.3600 s 0.4800 -0.1800 0.2000 0 -1.0400 1.2600 1.2800 0.5000 0.6800 -1.8600 0.6800 2.0000