三次样条插值,数值计算,上机报告,c.doc

三次样条插值,数值计算,上机报告,c 数值计算上机报告(三次样条) 1、题目 x0 1 4 9 16 25 36 49 64 y0 1 2 3 4 5 6 7 8 以上数据点的插值可以得到平方根函数的近似，在区间[0,64]上作图 用三次样条（第一边界条件）程序求S(x) 2、解题方法 在TC中建立一维数组，用以记录X, Y, h, μ, λ,d各值，建立二维数组用以记录三对角矩阵。将数据排列在三对角矩阵中，根据高斯消去法的原理解方程，算出各M值。通过解出的M值和h值计算三次样条的方程。 3、计算结果 区间[0,1], [1,4], [4,9], [9,16], [16,25], [25,36], [36,49], [49,64]上的三次样条函数分别是 8.915972*(1-x)*(1-x)*(1-x)-1.165277*x*x*x-8.915972*(1-x)+2.165277*x -0.388426*(4-x)*(4-x)*(4-x)+0.082175*(x-1)*(x-1)*(x-1)+3.829164*(4-x)-0.072911*(x-1) 0.049305*(9-x)*(9-x)*(9-x)-0.009943*(x-4)*(x-4)*(x-4)-0.832629*(9-x)+0.848581*(x-4) -0.007102*(16-x)*(16-x)*(16-x)+0.002111*(x-9)*(x-9)*(x-9)+0.776585*(16-x)+0.468011*(x-9) 0.001642*(25-x)*(25-x)*(25-x)-0.000168*(x-16)*(x-16)*(x-16)+0.311479*(25-x)+0.569198*(x-16) -0.000138*(36-x)*(36-x)*(36-x)+0.000154*(x-25)*(x-25)*(x-25)+0.471219*(36-x)+0.526880*(x-25) 0.00013*(49-x)*(49-x)*(49-x)+0.000074*(x-36)*(x-36)*(x-36)+0.439587*(49-x)+0.525918*(x-36) 0.000064*(64-x)*(64-x)*(64-x)-0.000041*(x-49)*(x-49)*(x-49)+0.452193*(64-x)+0.542654*(x-49) C程运行结果截图 Matlab画图结果截图 4、结果分析 运算过程中f #39; 的值不存在，当x 趋向于0 时， f x #39; 趋近于正无穷，因 此f #39; 取值不准确，在计算中， f #39; 取了100，因此前几个区间上的三 次样条误差较大，而随着自变量的增大，误差减小。与拉格朗日插值 相比三次样条插值对于此题的估计更为准确。如下图所示，绿色线条 表示的是八次拉格朗日插值，尽管在x 较小时拟合比较准确，但随着 自变量的增大，误差明显变大，稳定性较差。 5、源程序 #includestdio.h main() { double x[9]={0,1,4,9,16,25,36,49,64}; double y[9]={0,1,2,3,4,5,6,7,8}; double t, s, matrix[9][10], h[8], u[9], v[8], xishu[8][4]; int i, j; for(j=0;j=9;j++){ if(i==j) matrix[i][j]=2; else matrix[i][j]=0; } } for(i=0;i=7;i++){ h[i]=x[i+1]-x[i]; } u[0]=u[8]=1; for(i=1;i=7;i++){ u[i]=h[i-1]/(h[i-1]+h[i]); } v[0]=1; for(i=1;i=7;i++){ v[i]=h[i]/(h[i-1]+h[i]); } matrix[0][9]=100; for(i=1;i=7;i++){ matrix[i][9]=6*((y[i]-y[i+1])/(x[i]-x[i+1]))/(h[i-1]+h[i]); } matrix[i][9]=6*(1.0/16.0-(y[i-1]-y[i])/(x[i-1]-x[i]))/h[i-1]; for(j=0;j=8;j++){ for(i=0;i=8;i++){ if(i==(j+1)) matrix[j][i]=v[j]; else if(i==(j-1)) matrix[j][i]=u[j]; } } for(j=0;j=8;j++){ s=matrix[j][j]; t=matrix[j+1][j]; for(i=0;i=9;i++){ matrix[j][i